# ssfnong.ssc					script file for non-gaussian state space model support functions
#
# author: Eric Zivot
# created: September 24, 2002
# updates:
# September 25, 2002
#	Added comments and fixed a few typos
#	Added function SsfNonG.Antithetic - Siem Jan: please check for loop structure
# September 26, 2002
#	Added smaller functions SsfNonG.Expansion.SV, SsfNonG.LogDensity.SV
#	Eliminated loop from SsfNonG.Antithetic
# September 29, 2002
#	Added more comments to code
#	Fixed typos in SsfNonG.Expansion and SsfNonG.Expansion.SV
#
# notes:
# Code is adapted from OX code in ssfnong.ox
# available at www.ssfpack.com
#

SsfNonG.Expansion = function(iD,vY,vTheta=NULL,args) {
	##
	## Taylor series expansions for non-gaussian densities used
	## in state space models. See Durbin and Koopman (2000) JRSSB
	## inputs:
	## iD			character id for density. Valid choices are
	## 				d.normal, d.poisson, d.sv, d.svm, d.exp, d.t, 
	##				d.total
	## vY			vector of data
	## vTheta		vector of signals Z*alpha
	## args		vector optional arguments for densities
	## outputs:
	## mret		list with components y.tilde and H.tilde 
	##				from eqn (27) of Durbin and Koopman
	##
	vY = as.matrix(vY)
	vTheta = as.matrix(vTheta)
	iD = casefold(iD)		# convert to lower case characters
	if(iD == "d.normal") {
		y.tilde = vY
		H.tilde = matrix(1,nrow(vY),1)
	}
	else if(iD == "d.poisson") {
		if(missing(vTheta)) {
			y.tilde = log(vY)
			H.tilde = matrix(1,nrow(vY),1)
		}
		else {
			H.tilde = exp(-vTheta)
			y.tilde = vTheta + (H.tilde*vY) - 1
		}
	}
	else if(iD == "d.sv") {
		vy2 = vY*vY
		dsigma2 = args	# average variance
		if(missing(vTheta)) {
			y.tilde = log(vy2)-log(dsigma2)
			H.tilde = matrix(1,nrow(vY),1)
		}
		else {
			H.tilde = 2*dsigma2*exp(vTheta)/vy2	
			y.tilde = vTheta + 1 - 0.5*H.tilde
		}
	}
	else if(iD == "d.svm") {
		vy2 = vY*vY
		dsigma2 = args[1]
		dbeta2 = sqrt(args[2])
		if(missing(vTheta)) {
			y.tilde = log(vy2)-log(dsigma2-dbeta2)
			H.tilde = matrix(1,nrow(vY),1)
		}
		else {
			vx = exp(vTheta)
			vx2 = vx*vx
			vdp = vy2 + dbeta2*vx2
			vdm = vy2 - dbeta2*vx2
			H.tilde = (2*dsigma2*vx)/vdp
			y.tilde = vTheta + (vdm - dsigma2*vx)/vdp
		}
	}
	else if(iD == "d.exp") {
		if(missing(vTheta)) {
			y.tilde = log(vY)
			H.tilde = matrix(1,nrow(vY),1)
		}
		else {
			H.tilde = exp(vTheta)/vY
			y.tilde = vTheta + 1 - H.tilde
		}
	}
	else if (iD == "d.t") {
		dvar = args[1]
		sd = args[2]
		if(missing(vTheta)) {
			y.tilde = vY
			H.tilde = dvar
		}
		else {
			y.tilde = vY
			H.tilde = dvar*(df - 2) + (vTheta*vTheta)/(df + 1)
		}
	}
	else stop("Invalid choice for iD")
	mret = list(y.tilde = y.tilde, H.tilde = H.tilde)
	mret
}

SsfNonG.Expansion.SV = function(iD,vY,vTheta = NULL,args) {
	##
	## Taylor series expansions for non-gaussian densities used
	## in SV state space model. See Durbin and Koopman ch. 11, pg. 195
	## inputs:
	## iD			character id for density. Valid choices are
	## 				d.normal, d.sv
	## vY			vector of data
	## vTheta		vector of signals Z*alpha
	## args		vector optional arguments for densities
	## outputs:
	## mret		list with components y.tilde and H.tilde from 
	##				pg. 195 of Durbin and Koopman (2001)
	##
	vY = as.matrix(vY)
	vTheta = as.matrix(vTheta)
	iD = casefold(iD)	
	if(iD == "d.normal") {
		y.tilde = vY
		H.tilde = matrix(1,nrow(vY),1)
	}
	else if(iD == "d.sv") {
		vy2 = vY*vY
		dsigma2 = args	# average variance for SV model
		if(missing(vTheta)) {
			y.tilde = log(vy2) - log(dsigma2)
			H.tilde = matrix(1,nrow(vY),1)
		}
		else {
			H.tilde = 2*dsigma2*exp(vTheta)/vy2
			y.tilde = vTheta + 1 - 0.5*H.tilde
		}
	}
	else stop("Invalid choice for iD")
	mret = list(y.tilde = y.tilde, H.tilde = H.tilde)
	mret
}

SsfNonG.LogDensity = function(iD,mY,vTheta,args) {
	##
	## create Gaussian and non-Gaussian log densities for state space models
	## Gausian density is g(y|theta) and non-Gaussian density is p(y|theta)
	## These values are use for computing importance sampling weights
	##
	## inputs:
	## iD			character id for density. Valid choices are
	## 				d.normal, d.poisson, d.sv, d.svm, d.exp, d.t, 
	##				d.total
	## mY			vector of data
	## vTheta		vector of signals Z*alpha
	## args		vector optional arguments for densities
	## outputs:
	## dret		log-density lnp(y|theta) or lng(y|theta) without 
	##				normalizing constant term
	##
	mY = as.matrix(mY)
	vTheta = as.matrix(vTheta)
	iD = casefold(iD)
	if(iD == "d.normal") {
		dvar = args	# can be a vector for heteroskedastic model
		dret = -0.5*sum(log(dvar) + (vTheta*vTheta/dvar))
	}
	else if(iD == "d.poisson") {
		dret = crossprod(vTheta,mY) - sum(exp(vTheta))
	}
	else if(iD == "d.svm") {
		dsigma2 = args[1]
		vy = mY - (args[2] + (args[1]*exp(vTheta)) + args[3])
		dret = -0.5*(nrow(mY)*log(dsigma2) + sum(vTheta) +
				sum(vy*vy*exp(-vTheta))/dsigma2)
	}
	else if(iD == "d.sv") {
		dsigma2 = args
		dret = -0.5*(nrow(mY)*log(dsigma2) + sum(vTheta) +
				sum(mY*mY*exp(-vTheta))/dsigma2)
	}
	else if(iD == "d.exp") {
		dret = -sum(vTheta + mY*exp(-vTheta))
	}
	else if(iD == "d.t") {
		dvar = args[1]
		df = args[2]
		dk = (df - 2)*dvar
		dret = nrow(mY)*(lgamma((df + 1)/2) - lgamma(df/2) -
		0.5*log(dk)) - 0.5*(df + 1)*sum(log(1 + ((vTheta*vTheta)/dk)))
	}
	else stop("Invalid choice for iD")
	dret		# log-likelihood without constant term
}

SsfNonG.LogDensity.SV = function(iD,mY,vTheta,args) {
	##
	## create non-Gaussian log densities for SV state space models
	## non-Gaussian density is p(y|theta); Gaussian density is g(y|theta)
	## These values are used for computing importance sampling weights
	##
	## See Durbin and Koopman (2000) pg. 195
	## inputs:
	## iD			character id for density. Valid choices are
	## 				d.normal, d.sv
	## mY			vector of data
	## vTheta		vector of signals Z*alpha when id="d.sv"; vector of
	##				errors when id="d.normal"
	## args		scalar or vector of optional arguments for densities
	## outputs:
	## dret		log-density lnp(y|theta) or lng(y|theta) without 
	##				normalizing constants
	##
	mY = as.matrix(mY)
	vTheta = as.matrix(vTheta)
	iD = casefold(iD)
	if(iD == "d.normal") {
		dvar = args	# can be a vector for heteroskedastic model
		dret = -0.5*sum(log(dvar) + (vTheta*vTheta/dvar))
	}
	else if(iD == "d.sv") {
		dsigma2 = args	# avg variance
		dret = -0.5*(nrow(mY)*log(dsigma2) + sum(vTheta) +
				sum(mY*mY*exp(-vTheta))/dsigma2)
	}
	else stop("Invalid choice for iD")
	dret		# log-likelihood without constant term
}

SsfNonG.Antithetic = function(nrAnti, vDraw, vHat, dChi2) {
	##
	## Monte Carlo integration using anti-thetic variates
	##
	## inputs:
	## nrAnti			number of anti-thetic variates: 1,2 or 4
	## vDraw			vector of random draws from importance density
	## vHat			vector of smoothed disturbances
	## dChi2			degrees of freedom for chi-square
	## output:
	## mret			matrix of anti-thetic variates with nrAnti cols
	##
	## notes: 			See Durbin and Koopman JRSSB(2000) pg. 13 or
	##					Durbin and Koopma (2001), chs. 11 and 12
	if(!nrAnti==1 & !nrAnti==2 & !nrAnti==4)
		stop("nrAnti must be 1, 2 or 4")
	vHat = as.matrix(vHat)
	vDraw = as.matrix(vDraw)
	n = nrow(vHat)
	err = vDraw
	if (nrAnti == 1) {
		mret = err
	}
	else if(nrAnti == 2) {
		mret = cbind(err,2*vHat - err)
	}	
	else {
		d = sqrt(qchisq(1-pchisq(dChi2,n),n)/dChi2)
		err1 = vHat + d*(err - vHat)
		err2 = vHat + d*(vHat - err)
		mret = cbind(err,2*vHat - err,err1,err2)
	}
	mret
}