# fit SV model to NYSE returns  - uses SVfilter

y=matrix(scan("/mydata/nyse.dat"),ncol=1)
num=length(y)
y=log(y^2)

#initial parameters
phi0= 0
phi1=.95
sQ=.2
alpha=mean(y)  
sR0=1
mu1=-3
sR1=2
 initpar=c(phi0,phi1,sQ,alpha,sR0,mu1,sR1)

#-- evaluate the innovations likelihood
Linn=function(para){
phi0=para[1]
phi1=para[2]   
sQ=para[3]
alpha=para[4]
sR0=para[5]
mu1=para[6]   
sR1=para[7]
sv = SVfilter(num,y,phi0,phi1,sQ,alpha,sR0,mu1,sR1)
return(sv$like)
}

#-- do the estimation 
 est=optim(initpar, Linn, NULL, method = "BFGS", hessian = TRUE, control=list(trace=1,REPORT=1))
 stderr=sqrt(diag(solve(est$hessian)))
 est   # for a summary
 cbind(est$par,stderr)  # list estimates and SEs


#-- a nice graph --
phi0=est$par[1]
phi1=est$par[2]   
sQ=est$par[3]
alpha=est$par[4]
sR0=est$par[5]
mu1=est$par[6]   
sR1=est$par[7]
sv = SVfilter(num,y,phi0,phi1,sQ,alpha,sR0,mu1,sR1)

time=801:1000    # restrict time to 200 pts for viewing pleasure
par(mfrow=c(2,1))
plot(time, y[time], type="l", main="log(Squared NYSE Returns)", ylab="")
plot(time, sv$xp[time], type="l", main="Predicted log-Volatility", ylim=c(-1.5,1.8), ylab="")
lines(time,sv$xp[time]+2*sqrt(sv$Pp[time]), lty="dashed")
lines(time, sv$xp[time]-2*sqrt(sv$Pp[time]), lty="dashed")