################################################### # -- must source ss2 first -- ################################################### # -- read data -- newdat = read.table("/mydata/Newbold.dat", comment.char="#") y = newdat[1:50,2] # quarterly inflation (column 2) z = newdat[1:50,3] # quarterly interest rates (column 3) # # -- set up ----------------------- num=length(y) A=array(z, dim=c(1,1,num)) input=matrix(1,50,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 # nboot=200 # number of bootstrap replicates - #--------------------------------- # -- 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) } #-- # -- 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] 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=c(.84,-.77,.85,.12,1.1) # intial parameters for bootstrap (same as in estimation) for (i in 1:nboot){ cat("iteration:", i, "\n") e.star[4:50] = sample(e[4:50], replace=TRUE) # resample standardized innovations starting at t=4 for (j in 4:50){ xp.star[j] = phi*xp.star[j-1] + Ups + K[j]*sqrt(sig[j])*e.star[j] # generate bootstrapped states t=4,..,50 } y.star[4:50]=z[4:50]*xp.star[4:50] + alpha + sqrt(sig[4:50])*e.star[4:50] # 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- # reltol=.001 works well if you have the time # est.star=optim(initpar,Linn2,NULL,method="BFGS", control=list(reltol=.1)) # 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 # ----------------------------------------------------------------------------------- #-- 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,mux,Sigmax,phi,Ups,alpha,cQ,cR,0,input) # plot.ts(y) # lines(z, 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) # ---------------------------------------------------------------------------------------