## checkSV.ssc				script files to check SV estimation functions
##
## author: Eric Zivot
## created: September 29, 2002
## updates:
## May 26, 2003
##		Fixed problem with computing smoothed estimates of volatility from MSL estimation of SV model for 
##		exchange rates
## May 18, 2003
##		Added SML estimation of exchange rates. Results are close to but not equal to those available at
##		http://staff.feweb.vu.nl/koopman/sv/
##
##	May 12, 2003
##		Added time plots
##	
## notes:
##	Uses functions from ssfNong.ssc and sv_mcl_est.ssc
##	Data are the pound/Dollar daily exchange rates from 1/10/81 to 28/6/85. There are 945 obvs.
##	Raw data on log returns are in the data.frame svpdx.dat
##

## plot data UK/UK daily exchange rate data
## compare with plots on Siem-Jan's webpage
##
tsplot(svpdx.dat,main="Log returns on UK/US spot exchange rate",ylab="log return * 100",
xlab="day")
tsplot(abs(svpdx.dat),main="Absolute log returns of US/UK spot rate",ylab="Abs(log return*100)",
xlab="day")

## use parameter estimated from Durbin and Koopman (2001) page 236
##
s.dScale = (0.6338)^2
sig2.n = (0.1726)^2
phi = 0.9731

## use parameters from Siem-Jan's webpage
##
## s.dScale = 0.403462
## sig2.n = 0.027807
## phi = 0.9744


###########################################################################################################
## check out various functions
###########################################################################################################

## create state space form for linear gaussian approximating model
## AR(1) model with heteroskedastic measurement equation errors
## error variances will be computed later
##
ssf.lgam = GetSsfArma(ar=phi,sigma=sqrt(sig2.n))
s.mIndex = matrix(-1,2,2)
s.mIndex[2,2] = 1	# index for measurement eqn error

## starting values for iteration to find linear Gaussian
## approximating model (LGAM)
##
s.mY = as.matrix(svpdx.dat)
s.vTheta = matrix(0,nrow(s.mY),1)

## compute starting values for y.tilde and H.tilde for LGAM
## based on 2nd order Taylor series expansion
##
mret = SsfNonG.Expansion.SV(iD="d.sv",vTheta=s.vTheta,vY=s.mY,args=s.dScale)

## adjust state space to have time varying variances
##
ssf.lgam$mJOmega = s.mIndex
ssf.lgam$mX = mret$H.tilde

## compute approximating model
##
mApprox = SsfApproxModel("d.sv",s.mY,s.vTheta,s.dScale,ssf.lgam,s.mIndex,
					n.iter=50,tol=10^-7)


## evaluate log-likelihood using MC integration
##
set.seed(123)
loglik = SsfMCLik(s.mY,mApprox,s.dScale,ssf.lgam,s.mIndex,
s.iM=10,s.iAnti=1)
loglik

start.seed = 123
loglik2 = loglike.fun.SV(parm.start,mY=s.mY,start.seed)
loglik2

## estimate mean and variance of log-volatility
## using Monte Carlo Integration
ans = DrawMean(s.mY=s.mY,s.vTheta=s.vTheta,s.dScale=s.dScale,ssf=ssf.lgam,s.mIndex=s.mIndex,
				n.iter=50,tol=10^-7,	s.iM=250,s.iSeed=123) 

## plot smoothed estimates of volatility 
##
tsplot(ans$sd.hat,abs(svpdx.dat),ans$lower,ans$upper,
lty=c(1,2,1,1),col=c(3,1,2,2))
legend(0,3,legend=c("volatility","abs(returns)","-2SD","+2SD"),
lty=c(1,2,1,1),col=c(3,1,2,2))

ssf.lgam$mX = mApprox$H.tilde
kf.lgam = KalmanFil(mY=mApprox$y.tilde,ssf=ssf.lgam,task="STSIM")
m.s = SimSmoDraw(kf=kf.lgam,ssf=ssf.lgam,task="STSIM")

#############################################################################################################
## estimate SV model by SML estimation
#############################################################################################################
p1 = log(phi/(1-phi))
p2 = log(sqrt(sig2.n))
p3 = log(sqrt(s.dScale))
parm.start = c(p1,p2,p3)

## estimate SV model by SML
## using nlminb. results are close but not quite equal to those on Siem-Jan's webpage
##
tmp = nlminb.trace(parm.start,loglike.fun.SV,mY=s.mY,s.seed=start.seed,trace=T)
tmp.mle = tmp$parameters
phi.hat = exp(tmp.mle[1])/(1+exp(tmp.mle[1]))
sig2.n.hat = exp(2*tmp.mle[2])
s.dScale.hat = exp(2*tmp.mle[3])

ssf.lgam.hat = GetSsfArma(ar=phi.hat,sigma=sqrt(sig2.n.hat))
s.mIndex = matrix(-1,2,2)
s.mIndex[2,2] = 1	# index for measurement eqn error

## estimate mean and variance of log-volatility
## using Monte Carlo Integration
ans.hat = DrawMean(s.mY=s.mY,s.vTheta=s.vTheta,s.dScale=s.dScale.hat,ssf=ssf.lgam.hat,s.mIndex=s.mIndex,
				n.iter=50,tol=10^-7,	s.iM=250,s.iSeed=123) 

## plot smoothed estimates
## Note: picture does not look as nice as the one on Siem-Jan's web-page
##
tsplot(ans.hat$sd.hat,abs(svpdx.dat),ans.hat$lower,ans.hat$upper,
lty=c(1,2,1,1),col=c(3,1,2,2))
legend(0,3,legend=c("volatility","abs(returns)","-2SD","+2SD"),
lty=c(1,2,1,1),col=c(3,1,2,2))



##
## analysis of simulated SV data
##

## function to simulate from simple SV model
## y(t) = s*exp(0.5*theta(t))*u(t), u(t) ~ N(0,1)
## theta(t) = phi*theta(t-1) + e(t), e(t) ~ N(0,s2.e)
##
SV.sim = function(nsim=100,s.ave=1,phi=0.9,s.e=0.1) {
	## nsim 		number of simulations
	## s.ave		average sd
	## phi			autoregressive parameter
	## s.e			log volatility innov sd
	u = rnorm(nsim)
	e = rnorm(nsim,sd=s.e)
	theta = arima.sim(nsim,model=list(ar=phi),innov=e)
	y = s.ave*exp(0.5*theta)*u
	ans = list(y=y,theta=theta)
	ans
}

## simulate data calibrated to UK/US data
##
set.seed(333)
s.ave = 0.6338
phi = 0.9731
s.e = 0.1726
sim.SVdata = SV.sim(nsim=1000,s.ave=s.ave,phi=phi,s.e=s.e)
tsplot(sim.SVdata$y)
tsplot(abs(sim.SVdata$y),s.ave*exp(0.5*sim.SVdata$theta))