###################################################
#  -- must source ss2 first --
source("http://www.stat.pitt.edu/stoffer/tsa2/Rcode/New6/ssall.R")
###################################################
# -- read data --
newdat = read.table("http://www.stat.ualberta.ca/~wiens/stat479/S&Sdatasets/Newbold.dat", comment.char="@")
#z = newdat[1:50,2]   # quarterly inflation (column 2)
#y = newdat[1:50,3]   # quarterly interest rates (column 3)

z = newdat[1:100,2]   # quarterly inflation (column 2)
y = newdat[1:100,3]   # quarterly interest rates (column 3)

y1 = ts(y, start = 1953, frequency=4)
z1 = ts(z, start = 1953, frequency=4)
ts.plot(y1,z1,lty=1:2)  # Interest rate (y) is solid, inflation rate (z) is dashed; 
points(z1, pch=22)
points(y1, pch=16)

#
#-- number of bootstrap replicates - 
 nboot=200 ## Set this to a much smaller value at first, to test the algorithm         
#
#                     
# -- Function to Calculate Likelihood --
 Linn=function(para){
  phi=para[1]
  alpha=para[2]
  b=para[3]     
     Ups=(1-phi)*b
  cQ=para[4]      
  cR=para[5]  
  kf=Kfilter2(num,y,A,mu1,Sigma1,phi,Ups,alpha,cQ,cR,0,input)
  return(kf$like)
  }
#  
#
# -- Function to Calculate Likelihood for the bootstrapped data--
# -- Linn and Linn2 are the same except one uses y (the data)
#         and one uses y.star (the bootstrapped data)
 Linn2=function(para){
  phi=para[1]
  alpha=para[2]
  b=para[3]     
     Ups=(1-phi)*b
  cQ=para[4]      
  cR=para[5]  
  kf=Kfilter2(num,y.star,A,mu1,Sigma1,phi,Ups,alpha,cQ,cR,0,input)
  return(kf$like)
  }
  
# -- set up --
 num=length(y) 
 A=array(z, dim=c(1,1,num))
 input=matrix(1,num,1)     
 mu1=1; Sigma1=.01;             
 phi=.84; alpha=-.77;  b=.85; cQ=.12; cR=1.1
 initpar=c(phi,alpha,b,cQ,cR)    # initial parameters 

  
  
# --  bootstrap procedure --
  est=optim(initpar,Linn,NULL,method="BFGS", control=list(reltol=.001)) # estimate parameters
   phi=est$par[1];  alpha=est$par[2]; b=est$par[3]; Ups=(1-phi)*b   # the estimates
   cQ=est$par[4]; cR=est$par[5] 


# ----------------------------------------------------------------------------------- 
#  I copied this from the bottom of this script; 
#  it gives the same MLEs but gves the standard errors too
#------------------------------------------------------------------------------------
 est=optim(initpar,Linn,NULL,method="BFGS",hessian=TRUE,control=list(trace=1,REPORT=1))
 stderr=sqrt(diag(solve(est$hessian)))
 est      # to display summary of estimation 
   estimate=est$par
   u=cbind(estimate,stderr)
   rownames(u)=c("phi","alpha","b","sigw","sigv")
   u      # display parameter estimates and SEs

#------------------------------------------------------------------------------------

  kf=Kfilter2(num,y,A,mu1,Sigma1,phi,Ups,alpha,cQ,cR,0,input)       # run the filter
# pull out necessary values from the filter
   xp=kf$xp
   innov=kf$innov
   sig=kf$sig
   K=kf$K
   e=innov/sqrt(sig)
   e.star=e   # initialize resampled errors (just so e.star has the proper dimensions)
   y.star=y   # initialize data  -- the first 3 y values will be held fixed
   xp.star=xp # initialize state predictions
   para.star = matrix(0, nboot, 5)  # a place to put the bootstrapped estimates
   initpar = u[,1]# Start with the MLEs 
   for (i in 1:nboot){                             
      cat("iteration:", i, "\n")     
      e.star[4:num] = sample(e[4:num], replace=TRUE)   # resample standardized innovations starting at t=4
      for (j in 4:num){
       xp.star[j] = phi*xp.star[j-1] + Ups + K[j]*sqrt(sig[j])*e.star[j]  # generate bootstrapped states t=4,..,num
      } 
      y.star[4:num]=z[4:num]*xp.star[4:num] + alpha +  sqrt(sig[4:num])*e.star[4:num]  # generate bootstrapped obs t=4,..,50
# estimation with booted data... note reltol=.1 gets you stinky estimates but saves time-
#   if you set reltol to a lower value, you may be waiting awhile to complete this analysis      
   est.star=optim(initpar,Linn2,NULL,method="BFGS", control=list(reltol=.01))  # reltol set here
   para.star[i,]=cbind(est.star$par[1], est.star$par[2], est.star$par[3],   # record estimates
                 est.star$par[4], est.star$par[5])                      
   } 
# 
# the results are in para.star... phi, alpha, b, sigw, sigv are the 5 columns
# so for example
hist(para.star[,1], 15, main="phi")  
# gives a histogram of the phi estimates
 
#Display the bootstrap estimates and standard errors
boot.se=vector()
boot.est = vector()
for(i in 1:5) {
	boot.est[i] = mean(para.star[,i])
	boot.se[i] = sqrt(var(para.star[,i]))
	}
out = cbind(boot.est, boot.se)
rownames(out)=c("phi","alpha","b","sigw","sigv")
round(out,3)

 
# Bootstrap ci on phi:
phistar = para.star[,1]
quantile(phistar,c(.05,.95)) 

# Bootstrap ci on sigmaw:
sigmawstar = para.star[,4]
quantile(sigmawstar,c(.05,.95)) 

win.graph()
hist(para.star[,4], 15, main="sigmaw")  
# gives a histogram of the sigmaw estimates

# ----------------------------------------------------------------------------------- 
#-- Full estimation of the data --
#------------------------------------------------------------------------------------
# est=optim(initpar,Linn,NULL,method="BFGS",hessian=TRUE,control=list(trace=1,REPORT=1))
# stderr=sqrt(diag(solve(est$hessian)))
# est      # to display summary of estimation 
#   estimate=est$par
#   u=cbind(estimate,stderr)
#   rownames(u)=c("phi","alpha","b","sigw","sigv")
#   u      # display parameter estimates and SEs
 
#-- a plot of the results --
#  phi=est$par[1];  alpha=est$par[2]; b=est$par[3]; Ups=(1-phi)*b
#  cQ=est$par[4]; cR=est$par[5] 
#  kf=Kfilter2(num,y,A,mu1,Sigma1,phi,Ups,alpha,cQ,cR,0,input)
#  plot.ts(y1, lty=1, ylim=c(-5,5))
#  lines(z1, lty=2)
#  lines(unlist(kf$xp), col=2)
#  lines(unlist(kf$xp)+2*sqrt(unlist(kf$Pp)), col=4)
#  lines(unlist(kf$xp)-2*sqrt(unlist(kf$Pp)), col=4)
# ---------------------------------------------------------------------------------------