# SNPexamples.ssc				Examples for SNP chapter
# 
# authors: Eric Zivot, Jiahui Wang, Ying Gu
# created: July 1, 2003
# revised: March 30, 2005
#
# comments
# 1. Example data
#		dmRet.dat 	exchange rate data
#		tbill.dat		weekly interest rate data (US 3 month, 1 yr and 10 yr)
#		ckls.ts		monthly interest rate data (US 3 month)
#		msft.ts		daily returns on microsoft
#		sp500.ts		daily returns on sp500
# 2. The SNP model orders are represented by the vector of tuning parameters
#		(L_u,L_g,L_r,L_p,K_z,I_z,K_x,I_x) = 
#		(ar,garch,arch,lagP,zPoly,trimZ,xPoly,trimX)

#
# Exchange rate data
#

# summarize and plot spot rate and forward discount data
start(dmRet.dat)
end(dmRet.dat)
colIds(dmRet.dat)
par(mfrow=c(2,1))
	plot(dmRet.dat[,"exch"],main="US/DM spot exchange rate",
	reference.grid=F)
	abline(h=0)
	plot(dmRet.dat[,"forward"],main="US/DM forward discount",
	reference.grid=F)
	abline(h=0)
par(mfrow=c(1,1))
summaryStats(dmRet.dat[,c("exch","forward")])

# plot sample ACF for spot and forward discount rate data
tmp = acf(dmRet.dat[,c("exch","forward")])

#
# Gaussian VAR models
#

# fit Gaussian AR(1) (10010000) to exchange rate data
# setting aside first 14 observations
fit.10010000 = SNP(dmRet.dat[,"exch"],model=SNP.model(ar=1),
                   n.drop=14)
class(fit.10010000)
fit.10010000
summary(fit.10010000)

tmp = SNP(dmRet.dat[,"exch"],model=SNP.model(ar=1),
n.drop=14,control=SNP.control(version=8.7))


# fit without location-scale transform
fit.ar1 = SNP(dmRet.dat[,"exch"],model=SNP.model(ar=1),
n.drop=14,LStransform=list(mu=0,cstat=1,pstat=1))
coef(fit.ar1)
# check with FinMetrics function OLS
fit.OLS = OLS(exch~tslag(exch),data=dmRet.dat[14:834,])
coef(fit.OLS)

# fit Gaussian VAR(1) (10010000) to exchange rate data
# setting aside first 14 observations
bfit.10010000 = SNP(dmRet.dat[,1:2],model=SNP.model(ar=1),
                    n.drop=14)
class(bfit.10010000)
bfit.10010000
summary(bfit.10010000)

# fit without location-scale transformation - traditional VAR
fit.var1 = SNP(dmRet.dat[,1:2],model=SNP.model(ar=1),
n.drop=14,LStransform=list(mu=c(0,0),cstat=diag(2),pstat=diag(2)))
coef(fit.var1)
# check with FinMetrics function VAR
# note: arrangement of coefficients from SNP is different from the 
# arrangement of the coefficients from VAR
fit.VAR = VAR(cbind(exch,forward)~ar(1),data=dmRet.dat[14:834,])
print(coef(fit.VAR), digits=2)


#
# semi-parametric VAR models
#

# fit semi-parametric AR(1) (10014000) to exchange rate data
# setting aside first 14 observations
fit.10014000 = SNP(dmRet.dat[,"exch"],model=SNP.model(ar=1,zPoly=4),
n.drop=14)
fit.10014000
summary(fit.10014000)

# note: if zPoly < 3, then the polynomial terms are set to zero!
fit.10011000 = SNP(dmRet.dat[,"exch"],model=SNP.model(ar=1,zPoly=1),
n.drop=14)
fit.10011000

# fit restricted model with coefficients on z1 and z3 equal to zero
fit.10014000$coef
fit.10014000$est
coef.start = fit.10014000$coef
est.start = fit.10014000$est
coef.start[c(2,4)] = 0
est.start[c(2,4)] = F
fitr.10014000 = SNP(dmRet.dat[,"exch"],model=SNP.model(ar=1,zPoly=4),
n.drop=14,coef=coef.start,est=est.start)
fitr.10014000


# plot SNP density along with normal
SNP.density(fit.10014000)

# fit model with logistic transformation
fit.10014000.logit = SNP(dmRet.dat[,"exch"],model=SNP.model(ar=1,zPoly=4),
n.drop=14,control=SNP.control(xTransform="logistic"))
fit.10014000.logit

# fit semi-parametric VAR(1) (10014000) to exchange rate data
# setting aside first 14 observations
bfit.10014000 = SNP(dmRet.dat[,1:2],model=SNP.model(ar=1,zPoly=4),
n.drop=14)
bfit.10014000
summary(bfit.10014000)

SNP.density(bfit.10014000)

# utilize trimZ component
# exclude highest order interactions: trimZ=1
bfit.10014100 = SNP(dmRet.dat[,1:2],model=SNP.model(ar=1,zPoly=4,trimZ=1),
n.drop=14)
bfit.10014100
# show all terms in P(z)
bfit.10014100$idxZ

# filter out all high order interactions
bfit.10014300 = SNP(dmRet.dat[,1:2],model=SNP.model(ar=1,zPoly=4,trimZ=3),
n.drop=14)
bfit.10014300

#
# nonlinear nonparametric VAR
#

# fit nonlinear nonparametric AR(1) (10014000) to exchange rate data
# setting aside first 14 observations
# let P(z,x) only depend on first power of x
fit.10014010 = SNP(dmRet.dat[,"exch"],model=SNP.model(ar=1,zPoly=4,xPoly=1),
n.drop=14)
fit.10014010
summary(fit.10014010)

# fit model with P(z,x) depending on x(t-1) = (y(t-1),y(t-2))'
fit.10024010 = SNP(dmRet.dat[,"exch"],
model=SNP.model(ar=1,zPoly=4,xPoly=1,lagP=2),
n.drop=14)
fit.10024010

# fit nonlinear nonparametric VAR(1) (10014010) to exchange rate data
# setting aside first 14 observations
bfit.10014010 = SNP(dmRet.dat[,1:2],model=SNP.model(ar=1,zPoly=4,xPoly=1),
n.drop=14)
bfit.10014010

# exclude interactions 
bfit.10014321 = SNP(dmRet.dat[,1:2],
model=SNP.model(ar=1,zPoly=4,xPoly=2,trimZ=3,trimX=1),
n.drop=14)
bfit.10014321


#
# arch/garch specifications
#

# fit Gaussian AR(1) with GARCH(1,1) errors (11110000) to exchange rate data
# setting aside first 14 observations
fit.11110000 = SNP(dmRet.dat[,"exch"],model=SNP.model(ar=1,arch=1,garch=1),
n.drop=14)
class(fit.11110000)
fit.11110000
summary(fit.11110000)

# multivariate garch models
# darch		diagonal arch
# rgarch		restrict all garch coeffs to be the same
# dvec			diagonal vec model
# note: dvec model is not supported in Gallant and Tauchen Fortran code

# fit bivariate VAR(1) GARCH(1,1) with diagonal arch and restricted garch
bfit.11110000.r2 = SNP(dmRet.dat[,1:2],model=SNP.model(ar=1,arch=1,garch=1),
control=SNP.control(darch=T,rgarch=T,dvec=F),
n.drop=14)
bfit.11110000.r2

# fit bivariate VAR(1) GARCH(1,1) with diagonal arch 
bfit.11110000.r1 = SNP(dmRet.dat[,1:2],model=SNP.model(ar=1,arch=1,garch=1),
control=SNP.control(darch=T,rgarch=F,dvec=F),
n.drop=14)
bfit.11110000.r1

# fit bivariate VAR(1) GARCH(1,1)  
bfit.11110000 = SNP(dmRet.dat[,1:2],model=SNP.model(ar=1,arch=1,garch=1),
control=SNP.control(darch=F,rgarch=F,dvec=F),
n.drop=14)
bfit.11110000

# fit bivariate VAR(1) GARCH(1,1) DVEC model
bfit.11110000.dvec = SNP(dmRet.dat[,1:2],
model=SNP.model(ar=1,arch=1,garch=1),
control=SNP.control(dvec=T),
n.drop=14)
bfit.11110000.dvec

# fit bivariate VAR(1) GARCH(1,1) DVEC with restricted GARCH
# note: darch is redundant for DVEC model
bfit.11110000.dvec2 = SNP(dmRet.dat[,1:2],
model=SNP.model(ar=1,arch=1,garch=1),
control=SNP.control(darch=F,rgarch=T,dvec=T),
n.drop=14)
bfit.11110000.dvec2

# fit bivariate VAR(1) GARCH(2,2) DVEC model with restricted GARCH
bfit.12210000.dvec = SNP(dmRet.dat[,1:2],
model=SNP.model(ar=1,arch=2,garch=2),
control=SNP.control(dvec=T,rgarch=T),
n.drop=14)
bfit.12210000.dvec

#
# SNP estimation with random re-starts
#

# first fit an AR(1) model
fit.10010000 = SNP(dmRet.dat[,"exch"],model=SNP.model(ar=1),
n.drop=14)
fit.10010000



# first fit AR(1)-ARCH(3) with no re-starts
fit.10310000 = SNP(dmRet.dat[,"exch"],model=SNP.model(ar=1),
iArch=3,coef=fit.10010000$coef,n.drop=14)
fit.10310000


fit.10310000b = SNP(dmRet.dat[,"exch"],model=SNP.model(ar=1,arch=3),
n.drop=14)

# Re-estimate 10310000 model using 25 random restarts with perturbed 
# parameters. Note: fOld=0.1 perturbs the parameters of the start model
# by 0.1*U[0,1] and fNew=0.1 perturbs the additional parameters of the 
# fitted model by 0.1*U[0,1]

fit.10310000.restart = SNP(dmRet.dat[,"exch"],model=SNP.model(ar=1,arch=3),
n.drop=14,control=SNP.control(fOld=0.1,fNew=0.1,n.start=25),trace=T)

fOld = c(0, 1e-5, 0, -1e-5, 1e-5, 1e-5, 1e-4, 0, -1e-4, 1e-4, 1e-4, 
          1e-3, 0, -1e-3, 1e-3, 1e-3, 1e-5, 0, -1e-2, 1e-2, 1e-2, 
          1e-1, 0, -1e-1, 1e-1, 1e-1, 1e0, 0, -1e0, 1e0, 1e0) 
fNew = c(0, 0, 1e-5, 1e-5, -1e-5, 1e-5, 0, 1e-4, 1e-4, -1e-4, 1e-4, 
          0, 1e-3, 1e-3, -1e-3, 1e-3, 0, 1e-2, 1e-2, -1e-2, 1e-2, 
          0, 1e-1, 1e-1, -1e-1, 1e-1, 0, 1e0, 1e0, -1e0, 1e0) 
n.start = c(0,rep(25,30))
fit.10310000.mruns = SNP(dmRet.dat[,"exch"],model=SNP.model(ar=1,arch=3),  
            control = SNP.control(n.start = n.start, seed = 011667, 
            fOld = fOld, fNew = fNew), trace = T)

#
# The expand function
#
args(expand)

# expand Gaussian AR(1) to Gaussian AR(1)-ARCH(3)
fit.10310000 = expand(fit.10010000,arch=3)
class(fit.10310000)
fit.10310000

# expand with random re-starts using previously defined tweak constants
# There are 751 re-start estimations
# estimation takes about 20 seconds on 1.6 Ghz Pentium M
control.new = fit.10010000$control
control.new$fold = fOld
control.new$fnew = fNew
control.new$nstart = n.start
fit.10310000.mruns = expand(fit.10010000,arch=3,
control=control.new,trace=T)
fit.10310000.mruns

# compare AR(1) with AR(1)-ARCH(3) model selection criteria
fit.10010000$obj - fit.10310000.mruns$obj

#
# The SNP.auto function
#
args(SNP.auto)

# use SNP.auto on spot returns
fit.auto = SNP.auto(dmRet.dat[,"exch"],n.drop=14)

# use SNP.auto on spot returns with random re-starts
fit.auto.restart = SNP.auto(dmRet.dat[,"exch"],n.drop=14,
           control = SNP.control(n.start = n.start, seed = 011667, 
                     fOld = fOld, fNew = fNew))

#
# SNP model diagnostics
#

# extractor functions
# coef
# residuals
# fitted

# illustrate extractor functions using Gaussian AR(1)-GARCH(1,1)

coef(fit.11110000)
coef.fit.11110000 = coef(fit.11110000)
names(coef.fit.11110000)
# note: coef component contains unscrambled coefficients as well as coefficient names
coef.fit.11110000$coef

resid.11110000 = residuals(fit.11110000)
resid.11110000[1:3]
resid.11110000.std = residuals(fit.11110000,stardard=T)
sigma.t.11110000 = sigma.t.SNP(fit.11110000)

# plot conditional standard deviations along with moving average
# of log returns
sig.11110000 = sigma.t.SNP(fit.11110000)
par(mfrow=c(2,1))
plot(SMA(abs(dmRet.dat[15:834,"exch"])),
     main="Moving average of absolute log returns",
     reference.grid=F)
plot(sqrt(sig.11110000),main="conditional standard deviations from SNP fit",
     reference.grid=F)
par(mfrow=c(1,1))


# plot standardized residuals
plot(residuals(fit.11110000,stardard=T),reference.grid=F)

# illustrate plot method for univariate data
plot(fit.11110000)
plot(fit.11110000,ask=F,which.plots=7)


# illustrate plot method for multivariate data
plot(bfit.10010000)

plot(bfit.10010000,ask=F,which.plots=c(3,7))

# generate simulation from univariate SNP models
sim.10010000 = simulate(fit.10010000)

# generate simulation from semi-parametric AR(1) - GARCH(1,1)
fit.11114000 = SNP(dmRet.dat[,"exch"],
                   model=SNP.model(ar=1,arch=1,garch=1,zPoly=4),
                   n.drop=14)
sim.11114000 = simulate(fit.11114000)

par(mfrow=c(3,1))
	plot(dmRet.dat[,"exch"],reference.grid=F,
	main="US/DM spot returns")
	abline(h=0)
	plot(sim.10010000,reference.grid=F,
	main="Simulated returns from Gaussian AR(1) SNP model")
	abline(h=0)
	plot(sim.11114000,reference.grid=F,
	main="Simulated returns from semi-parametric AR(1)-GARCH(1,1) SNP model")
	abline(h=0)
par(mfrow=c(1,1))


par(mfrow=c(2,1))
plot(dmRet.dat[,"exch"],reference.grid=F,
main="US/DM spot returns")
abline(h=0)
plot(sim.11114000,reference.grid=F,
main="Simulated returns from semi-parametric AR(1)-GARCH(1,1) SNP model")
abline(h=0)
par(mfrow=c(1,1))

#
# predict method
#

args(predict.SNP)

# compute predictions from AR(1) model and compare to FinMetrics

predict.11110000 = predict(fit.11110000,n.predict=10)
class(predict.11110000)
predict.11110000
summary(predict.11110000)
# plot with 95% confidence bands
plot(predict.11110000,dmRet.dat[,"exch"],n.old=10,width=2)


# predictions using FinMetrics
fit.ols = OLS(exch~tslag(exch),data=dmRet.dat[14:834,])
predict.ols = predict(fit.ols,n.predict=10,olddata=dmRet.dat[,"exch"])


#
# data transformations
#

# logistic transform

logisticTransform <- function(x,sig.tr=2) {
	x.tr = (4*sig.tr)*exp(x/sig.tr)/(1+exp(x/sig.tr)) - 2*sig.tr
	x.tr
}

splineTransform <- function(x,sig.tr=4) {
	x.tr = ifelse((x<(-1*sig.tr)),.5*(x-sig.tr-log(1-x-sig.tr)),x)
	x.tr = ifelse((x>sig.tr),.5*(x+sig.tr+log(1-sig.tr+x)),x.tr)
	x.tr
}

# plot transformations
x = seq(from=-6,to=6,length=100)
x.logit = logisticTransform(x,sig.tr=2)
x.spline = splineTransform(x,sig.tr=2)
plot(x,x.spline,type="n",xlab="x",ylab="")
abline(a=0,b=1,lwd=2)
points(x,x.spline,type="l",lty=2,lwd=2)
points(x,x.logit,type="l",lty=3,lwd=2)
abline(v=c(-2,2))
legend(-6,4,legend=c("spline","logit"),lty=c(2,3),lwd=2)
# illustrate spline transformation
tb3mo = tbill.dat[,"tb3mo"]
fit.204100 = SNP(tbill.dat[,"tb3mo"],model=SNP.model(ar=2,arch=4),
n.drop=6)

# simulate from estimated AR(2)-ARCH(4) model
sim.204100 = simulate(fit.204100)

# fit with spline transformation
fit.204100.s = SNP(tbill.dat[,"tb3mo"],model=SNP.model(ar=2,arch=4),
n.drop=6,
control=SNP.control(xTransform="spline",inflection=4))

# simulate from estimated AR(2)-ARCH(4) model
sim.204100.s = simulate(fit.204100.s)

# plot simulations
par(mfrow=c(3,1))
plot(tb3mo,main="T-bill data",reference.grid=F)
plot(sim.204100,main="simulated data from fitted SNP 204100 model",
reference.grid=F)
plot(sim.204100.s,main="simulated data from SNP 204100 model with spline transform",
reference.grid=F)

par(mfrow=c(1,1))

par(mfrow=c(2,1))
plot(tb3mo,main="T-bill data",reference.grid=F)
plot(sim.204100.s,main="simulated data from SNP 204100 model with spline transform",
reference.grid=F)
par(mfrow=c(1,1))

############################################################################
# Detailied examples
############################################################################

#
# Fit SNP model to daily returns on Microsoft stock
#

start(msft.ts)
end(msft.ts)
nrow(msft.ts)

# plot data and compute summary statistics
par(mfrow=c(2,1))
	plot(msft.ts,reference.grid=F)
	abline(h=0)
	tmp=acf(msft.ts^2)
par(mfrow=c(1,1))

summaryStats(msft.ts)

msft.mat = as.matrix(seriesData(msft.ts))
par(mfrow=c(2,2))
	hist(msft.mat,xlab="returns",
	probability=T)
	boxplot(msft.mat,outchar=T,ylab="returns")
	plot(density(msft.mat),type="l",xlab="returns",
	ylab="density estimate")
	qqnorm(msft.mat,ylab="returns")
	qqline(msft.mat)
par(mfrow=c(1,1))

# model selection
# xTransform=logistic
# random re-starts after AR(1) fit


fOld <- c(0,1e-5, 0, -1e-5, 1e-5, 1e-5, 1e-4, 0, -1e-4, 1e-4, 1e-4,
          1e-3, 0, -1e-3, 1e-3, 1e-3, 1e-2, 0, -1e-2, 1e-2, 1e-2,
          1e-1, 0, -1e-1, 1e-1, 1e-1, 1e0, 0, -1e0, 1e0, 1e0)
fNew <- c(0,0, 1e-5, 1e-5, -1e-5, 1e-5, 0, 1e-4, 1e-4, -1e-4, 1e-4,
          0, 1e-3, 1e-3, -1e-3, 1e-3, 0, 1e-2, 1e-2, -1e-2, 1e-2,
          0, 1e-1, 1e-1, -1e-1, 1e-1, 0, 1e0, 1e0, -1e0, 1e0)
n.start <- c(0,rep(25,30))

# manual search for best model
fit.msft.10010000 <- SNP(msft.ts, model=SNP.model(ar=1), 
                     control = SNP.control(xTransform="logistic", 
                                           n.start=n.start, 
                                           fOld=fOld, fNew=fNew),
                     n.drop=14, trace=T)
fit.msft.11110000 <- expand(fit.msft.10010000, arch=1, garch=1, trace=T)
fit.msft.11114000 <- expand(fit.msft.11110000, zPoly=4)
fit.msft.11116000 <- expand(fit.msft.11114000, zPoly=2)
fit.msft.11116010 <- expand(fit.msft.11116000, xPoly=1)
fit.msft.11116020 <- expand(fit.msft.11116010, xPoly=1)
fit.msft.11216000 <- expand(fit.msft.11116000, arch=1)
fit.msft.12116000 <- expand(fit.msft.11116000, garch=1)
fit.msft.11118000 <- expand(fit.msft.11116000, zPoly=2)
## the best model we found is fit.msft.11116000  

# fit restricted SNP model imposing symmetry
# Set starting value for sigma=1
fit.msft.11116s00 <- SNP(data = msft.ts, 
                     model = SNP.model(ar=1, arch=1, garch=1, zPoly=6), 
                     control = SNP.control(xTransform="logistic", 
                                           n.start=n.start, 
                                           fOld=fOld, fNew=fNew), 
                     n.drop = 14, trace=T,
                     coef=c(1,0,0,0,0,0,0,0,0,1,0,0), 
                     est = c(0,0,1,0,1,0,1,1,1,1,1,1))

# automatic model selection
fit.msft.auto = SNP.auto(msft.ts,n.drop=14,
                control=SNP.control(xTransform="logistic",
                                    n.start=n.start, seed=011667,
                                    fOld=fOld,fNew=fNew))
# best fit is 01116000
summary(fit.msft.auto)

# diagnostics of best fitting model

# density plot
plot(fit.msft.11116000, ask=F, which.plot=c(2,3,7))
# simulated data
sim.msft.11116000 = simulate(fit.msft.11116000)
summaryStats(sim.msft.11116000)


## simulated plot vs the real data
par(mfrow=c(2,1))
	plot(data2.MSFT, main="stock returns of MSFT", reference.grid=F)
	plot(sim.msft.11116000, main="simulated returns of MSFT", 
	     reference.grid=F)
par(mfrow=c(1,1))

plot(sim.msft.11116000, reference.grid=F)
abline(h=0)

# compute distribution diagnostics
sim.msft.mat = as.matrix(seriesData(sim.msft.ts))
par(mfrow=c(2,2))
	hist(sim.msft.mat, xlab="returns",
	probability=T)
	boxplot(sim.msft.mat, outchar=T, ylab="returns")
	plot(density(sim.msft.mat), type="l", xlab="returns",
	ylab="density estimate")
	qqnorm(sim.msft.mat, ylab="returns")
	qqline(sim.msft.mat)
par(mfrow=c(1,1))

#
# Fit SNP model to daily returns on S&P 500 index
#

start(sp500.ts)
end(sp500.ts)
nrow(sp500.ts)

# plot data and compute summary statistics

par(mfrow=c(2,1))
	plot(sp500.ts,reference.grid=F)
	abline(h=0)
	tmp=acf(sp500.ts^2)
par(mfrow=c(1,1))

summaryStats(sp500.ts)

tmp = acf(abs(sp500.ts))

sp500.mat = as.matrix(seriesData(sp500.ts))
	par(mfrow=c(2,2))
	hist(sp500.mat,xlab="returns",
	probability=T)
	boxplot(sp500.mat,outchar=T,ylab="returns")
	plot(density(sp500.mat),type="l",xlab="returns",
	ylab="density estimate")
	qqnorm(sp500.mat,ylab="returns")
	qqline(sp500.mat)
par(mfrow=c(1,1))

# find best fitting model using manual model selection
fit.sp.10010000 <- SNP(data2.SP, model=SNP.model(ar=1), 
                   control = SNP.control(xTransform="logistic", 
                                         n.start=n.start, 
                                         fOld=fOld, fNew=fNew),
                   n.drop=14, trace=T)
fit.sp.11110000 <- expand(fit.sp.10010000, arch=1, garch=1)
fit.sp.11114000 <- expand(fit.sp.11110000, zPoly=4)
fit.sp.11116000 <- expand(fit.sp.11114000, zPoly=2)
fit.sp.11118000 <- expand(fit.sp.11116000, zPoly=2)
fit.sp.1111a000 <- expand(fit.sp.11118000, zPoly=2)

##
fit.sp.11118010 <- expand(fit.sp.11118000, xPoly=1)
fit.sp.11118020 <- expand(fit.sp.11118010, xPoly=1)

## the best model we found is fit.sp.11118000 

# find best fitting model using automatic model selection
fit.sp500.auto = SNP.auto(sp500.ts,n.drop=14,
                 control = SNP.control(xTransform="logistic",
                                       n.start = n.start, 
                                       seed = 011667, 
                                       fOld = fOld, fNew = fNew))
# best fitting model is a SNP 01118000 model
fit.sp500.auto

# fit GT's preferred models - 20b14000, where b = 11
fOld <- c(0,1e-5, 0, -1e-5, 1e-5, 1e-5, 1e-4, 0, -1e-4, 1e-4, 1e-4,
          1e-3, 0, -1e-3, 1e-3, 1e-3, 1e-2, 0, -1e-2, 1e-2, 1e-2,
          1e-1, 0, -1e-1, 1e-1, 1e-1, 1e0, 0, -1e0, 1e0, 1e0)
fNew <- c(0,0, 1e-5, 1e-5, -1e-5, 1e-5, 0, 1e-4, 1e-4, -1e-4, 1e-4,
          0, 1e-3, 1e-3, -1e-3, 1e-3, 0, 1e-2, 1e-2, -1e-2, 1e-2,
          0, 1e-1, 1e-1, -1e-1, 1e-1, 0, 1e0, 1e0, -1e0, 1e0)
n.start <- c(0,rep(25,30))
fit.sp500.20010000 <- SNP(sp500.ts, model=SNP.model(ar=2), 
                          control = SNP.control(xTransform="logistic", 
                                                n.start=n.start, 
                                                fOld=fOld, fNew=fNew),
                          n.drop=14, trace=F)
fit.sp500.20110000 = expand(fit.sp500.20010000, arch=1, trace=F)
fit.sp500.20210000 = expand(fit.sp500.20110000, arch=1)
fit.sp500.20310000 = expand(fit.sp500.20210000, arch=1)
fit.sp500.20410000 = expand(fit.sp500.20310000, arch=1)
fit.sp500.20510000 = expand(fit.sp500.20410000, arch=1)
fit.sp500.20610000 = expand(fit.sp500.20510000, arch=1)
fit.sp500.20710000 = expand(fit.sp500.20610000, arch=1)
fit.sp500.20810000 = expand(fit.sp500.20710000, arch=1)
fit.sp500.20910000 = expand(fit.sp500.20810000, arch=1)
fit.sp500.20a10000 = expand(fit.sp500.20910000, arch=1)
fit.sp500.20b10000 = expand(fit.sp500.20a10000, arch=1)

fit.sp500.20b14000 = expand(fit.sp500.20b10000, zPoly=4)
fit.sp500.20b14010 = expand(fit.sp500.20b14000, xPoly=1)

# fit simple 10014000 model
fit.sp500.10014000 <- SNP(sp500.ts, model=SNP.model(ar=1, zPoly=4), 
                          control = SNP.control(xTransform="logistic", 
                                                n.start=n.start, 
                                                fOld=fOld, fNew=fNew),
                          n.drop=14, trace=T)


# diagnostics for best fitting model
# density, std residuals, ACF of squared std residuals
# note: there is still some correlation in squared std residuals
plot(fit.sp.11118000, ask=F, which.plot=c(2,3,7))
# simulated data
sim.sp.11118000 = simulate(fit.sp.11118000)

## simulated plot vs the real data

par(mfrow=c(2,1))
	plot(sp500.ts,main="Daily returns on S&P 500", reference.grid=F)
	abline(h=0)
	plot(sim.sp.11118000, main="simulated returns on SNP model",
	     reference.grid=F)
	abline(h=0)
par(mfrow=c(1,1))

sim.sp.mat = as.matrix(seriesData(sim.sp.ts))
par(mfrow=c(2,2))
	hist(sim.sp.mat,xlab="returns",
	probability=T)
	boxplot(sim.sp.mat,outchar=T,ylab="returns")
	plot(density(sim.sp.mat),type="l",xlab="returns",
	ylab="density estimate")
	qqnorm(sim.sp.mat,ylab="returns")
	qqline(sim.sp.mat)
par(mfrow=c(1,1))


#
# Fit SNP model to weekly T-bill data
#

# data is in timeSeries object tbill.dat taken from Gallant and Tauchen
# weekly US 3 month T-bill, annualized rate times 100. Similar to data
# used by Andersen and Lund (1997)

colIds(tbill.dat)

tb3mo = tbill.dat[,"tb3mo"]
	par(mfrow=c(2,1))
	plot(tb3mo,reference.grid=F,main="")
	tmp=acf(tb3mo)
par(mfrow=c(1,1))
summaryStats(tb3mo)

# show distribution diagnostics for residuals from AR(1) model
resid.ar1 = residuals(OLS(tb3mo~tslag(tb3mo), data=tbill.dat, na.rm=T))
summaryStats(resid.ar1)

tb3mo.mat = as.matrix(seriesData(resid.ar1))
	par(mfrow=c(2,2))
	hist(tb3mo.mat,xlab="percent times 100",
	probability=T)
	boxplot(tb3mo.mat,outchar=T,ylab="percent times 100")
	plot(density(tb3mo.mat),type="l",xlab="percent times 100",
	ylab="density estimate")
	qqnorm(tb3mo.mat)
	qqline(tb3mo.mat)
par(mfrow=c(1,1))

#
# fit SNP models
# 

# 1. find best fitting model using SNP.auto
fit.tb3mo.auto <- SNP.auto(tb3mo, 
                  control = SNP.control(xTransform = "spline", 
                                        fNew=fNew, fOld= fOld, 
                                        n.start=n.start), 
                  arMax = 6, zPolyMax = 8, xPolyMax = 4, lagPMax=1)
fit.tb3mo.auto
fit.tb3mo.auto$obj

# initial SNP fit
fit.tb3mo.10010000 <- SNP(tb3mo, model=SNP.model(ar=1), 
                    control = SNP.control(xTransform="spline",
                                          fOld = fOld, fNew = fNew,
                                          n.start = n.start),
                    n.drop=14)

# 2. fit AR(1)-ARCH models

fit.tb3mo.10110000 <- expand(fit.tb3mo.10010000, arch=1)
fit.tb3mo.10210000 <- expand(fit.tb3mo.10110000, arch=1)
fit.tb3mo.10310000 <- expand(fit.tb3mo.10210000, arch=1)
fit.tb3mo.10410000 <- expand(fit.tb3mo.10310000, arch=1)
fit.tb3mo.10510000 <- expand(fit.tb3mo.10410000, arch=1)
fit.tb3mo.10610000 <- expand(fit.tb3mo.10510000, arch=1)
fit.tb3mo.10710000 <- expand(fit.tb3mo.10610000, arch=1)
fit.tb3mo.10810000 <- expand(fit.tb3mo.10710000, arch=1)
fit.tb3mo.10814000 <- expand(fit.tb3mo.10810000, zPoly = 4)
fit.tb3mo.10815000 <- expand(fit.tb3mo.10814000, zPoly = 1)
fit.tb3mo.10816000 <- expand(fit.tb3mo.10815000, zPoly = 1)
fit.tb3mo.10817000 <- expand(fit.tb3mo.10816000, zPoly = 1)
fit.tb3mo.10818000 <- expand(fit.tb3mo.10817000, zPoly = 1)
fit.tb3mo.10818010 <- expand(fit.tb3mo.10818000, xPoly = 1)
fit.tb3mo.10818020 <- expand(fit.tb3mo.10818010, xPoly = 1)

# 3. fit AR(1)-GARCH(1,1) models

fit.tb3mo.11110000 <- expand(fit.tb3mo.10010000, arch=1, garch=1, trace=T)
fit.tb3mo.11114000 <- expand(fit.tb3mo.11110000, zPoly=4)
fit.tb3mo.11115000 <- expand(fit.tb3mo.11114000, zPoly=1)
fit.tb3mo.11116000 <- expand(fit.tb3mo.11115000, zPoly=1)
fit.tb3mo.11117000 <- expand(fit.tb3mo.11116000, zPoly=1)
fit.tb3mo.11118000 <- expand(fit.tb3mo.11117000, zPoly=1)
fit.tb3mo.11118010 <- expand(fit.tb3mo.11118000,xPoly=1)
fit.tb3mo.11118020 <- expand(fit.tb3mo.11118010,xPoly=1)

# 4. fit G&T's preferred models in EMM user's Guide

fit.tb3mo.10414000 <- expand(fit.tb3mo.10410000, zPoly=4)
fit.tb3mo.10414010 <- expand(fit.tb3mo.10414000, xPoly=1)

#
# diagnostics of 11118000
#

# simulate from SNP model
tb3mo.11118000.sim = simulate(fit.tb3mo.11118000)
plot(tb3mo.11118000.sim, reference.grid=F)

par(mfrow=c(2,1))
	plot(fitted(fit.tb3mo.11118000), main="Fitted values",
	     reference.grid=F)
	plot(sigma.t.SNP(fit.tb3mo.11118000), main="Conditional SD values",
	     reference.grid=F)
par(mfrow=c(1,1))

## conditional density:

SNP.density(fit.tb3mo.11118000, scale=0.1)

## ACF of Std. Residuals and Squared Std. Residuals

plot(fit.tb3mo.11118000, ask=F, which.plots=5)
plot(fit.tb3mo.11118000, ask=F, which.plots=7)

#
# diagnostics of 10818000
#

# simulate from SNP model
tb3mo.10818000.sim = simulate(fit.tb3mo.10818000)
plot(tb3mo.10818000.sim, reference.grid=F)

par(mfrow=c(2,1))
	plot(fitted(fit.tb3mo.10818000), main="Fitted values",
	     reference.grid=F)
	plot(sigma.t.SNP(fit.tb3mo.10818000), main="Conditional SD values",
	     reference.grid=F)
par(mfrow=c(1,1))

## conditional density:

SNP.density(fit.tb3mo.10818000, scale=0.1)

## ACF of Std. Residuals and Squared Std. Residuals

plot(fit.tb3mo.10818000, ask=F, which.plots=5)
plot(fit.tb3mo.10818000, ask=F, which.plots=7)


#
# diagnostics of GT's preferred models 
#

tb3mo.10414000.sim = simulate(fit.tb3mo.10414000)
tb3mo.10414010.sim = simulate(fit.tb3mo.10414010)

par(mfrow=c(2,1))
plot(tb3mo.10414000.sim, main="simulated data from fitted 10414000 model",reference.grid=F)
plot(tb3mo.10414010.sim,main="simulated data from fitted 10414010 model",reference.grid=F)
par(mfrow=c(1,1))

## estimated conditional standard deiations

sigma.10414000=sigma.t.SNP(fit.tb3mo.10414000)
sigma.10414010=sigma.t.SNP(fit.tb3mo.10414010)

par(mfrow=c(2,1))
plot(sqrt(sigma.10414000), reference.grid=F)
plot(sqrt(sigma.10414010), reference.grid=F)
par(mfrow=c(1,1))

## conditional density

SNP.density(fit.tb3mo.10414000, scale=0.1)
SNP.density(fit.tb3mo.10414010, scale=0.1)

## Std. residuals:

plot(residuals(fit.tb3mo.10414000,stardard=T),reference.grid=F)
plot(fit.tb3mo.10414000,ask=F, which.plots=3)
plot(residuals(fit.tb3mo.10414010,stardard=T),reference.grid=F)
plot(fit.tb3mo.10414010, ask=F, which.plots=3)

## ACF of Std. Residuals and Squared Std. Residuals

plot(fit.tb3mo.10414000, ask=F, which.plots=5)
plot(fit.tb3mo.10414000, ask=F, which.plots=7)

#
# diagnostics of SNP.auto
#

# simulate from SNP model
tb3mo.auto.sim = simulate(fit.tb3mo.auto)
plot(tb3mo.auto.sim, reference.grid=F)

par(mfrow=c(2,1))
	plot(fitted(fit.tb3mo.auto), main="Fitted values",
	     reference.grid=F)
	plot(sigma.t.SNP(fit.tb3mo.auto), main="Conditional SD values",
	     reference.grid=F)
par(mfrow=c(1,1))

## conditional density:

SNP.density(fit.tb3mo.auto, scale=0.1)

## ACF of Std. Residuals and Squared Std. Residuals

plot(fit.tb3mo.auto, ask=F, which.plots=5)
plot(fit.tb3mo.auto, ask=F, which.plots=7)


#
# compare 4 models
#

par(mfrow=c(4,2))
plot(tb3mo.11118000.sim, main="11118000 model", reference.grid=F)
abline(h=0)
plot(sigma.t.SNP(fit.tb3mo.11118000), main="11118000 model", 
     reference.grid=F)
plot(tb3mo.10818000.sim, main="10818000 model",reference.grid=F)
abline(h=0)
plot(sigma.t.SNP(fit.tb3mo.10818000), main="10818000 model",
     reference.grid=F)
plot(tb3mo.auto.sim, main="51118000 model",reference.grid=F)
abline(h=0)
plot(sigma.t.SNP(fit.tb3mo.auto), main="51118000 model",
     reference.grid=F)
plot(tb3mo.10414010.sim,main="10414010 model",reference.grid=F)
abline(h=0)
plot(sigma.t.SNP(fit.tb3mo.10414010), main="10414010 model",
     reference.grid=F)

par(mfrow=c(1,1))

#
# Fit SNP model to monthly CRSP 30-day US T-bill used by CKLS
# note: use Mike Cliff's data ckls.ts
#

# data, taken from FinMetrics timeSeries rf.30day
# use CRSP 30 day T-bill data to match data sample used in CKLS
# and Andersen and Sorensen

# note: ckls.ts data is used in GMM chapter
par(mfrow=c(2,1))
	plot(ckls.ts, reference.grid=F, main="")
	tmp=acf(ckls.ts)
par(mfrow=c(1,1))
summaryStats(ckls.ts)

par(mfrow=c(2,1))
	plot(diff(ckls.ts), reference.grid=F)
	abline(h=0)
	tmp=acf(diff(ckls.ts))
par(mfrow=c(1,1))




# diagnostics of preferred model
# plot fitted values and conditional sd values (as in Tauchen (1997))
# Plot simulations from model


######################################################################################################
# dump SNP objects
######################################################################################################

dmp.names = c("fit.sp500.11118000","fit.sp500.auto","fit.sp500.20b14000",
              "fit.sp500.20b14010","fit.sp500.10014000","fit.msft.11116000",
              "fit.msft.auto")
data.dump(dmp.names,file="snpObjects.dmp")