# garch.ssc			
# author: Eric Zivot
# created: March 13, 2002
# updated: April 19, 2002

# Commented Version by JMS November 8, 2004.
# Saved as Ch7.garch.txt


module(finmetrics)

#Ford Daily Returns for 2000 days.

class(ford.s)
ford.s@title
start(ford.s)
end(ford.s)
nrow(ford.s)

#A look at the data

plot(ford.s,reference.grid=F)

#Now a look at the ACF and the ACF of the SQUARES.

par(mfrow=c(1,2))
tmp = acf(ford.s,lag=12)
tmp = acf(ford.s^2, lag=12)
par(mfrow=c(1,1)) 

#We see that the SQUARES have more lags of significance.
#To some extent this is a "cheat" since it is not clear 
#that the squared data will conform to the assumption that 
#one needs to set significance levels on the estimated correlations.
#Still the picture is suggestive.

#######################################
#Simulation of an ARCH(1) model

sim.arch1=simulate.garch(model=list(a.value=0.01, arch=0.8), n=250, rseed=196)
names(sim.arch1)

#The simulation function returns both the error term and the conditional standard 
#deviations. The ratio should be perfect normals.

ratio=sim.arch1$et/sim.arch1$sigma.t
normalTest(ratio, method="sw")
normalTest(ratio, method="jb")

#####################################
# Back to the real data and "testing for arch"
# More precisely we test that a_0=a_1=...=a_p where p=lag.n

archTest(ford.s, lag.n=12)

#We get a brutally small p-value, so it looks like there is an arch effect.
#Let estimate a univarite GARCH model for the series.
#What makes the model "Generalized" is that sigma.t is 
#modeled by lags in sigma.t as well as lags in epsilon.t
#We'll try one of each.

ford.mod11 = garch(ford.s~1, ~garch(1,1))
class(ford.mod11)
names(ford.mod11)

#Look at the default print...
ford.mod11

#More instructive to look at the summary
summary(ford.mod11)

#This has some weird inferences...
#Shapiro-Wilk is nonsignificant but Jarque-Berra is highly significant.
#This "may" mean that the residuals are "normal in the middle" but
#exhibit skewness or kurtosis. Also, the JB test is VERY sensitive
#to an outlier the series.

##############################################################
#Now lets look at some individual components of the series.

ford.mod11$asymp.sd	
	
coef(ford.mod11)
vcov(ford.mod11)
vcov(ford.mod11,method="qmle")
residuals(ford.mod11,standardize=T)
sigma.t(ford.mod11)

# garch diagnostics

summary(ford.mod11)
autocorTest(residuals(ford.mod11,standardize=T),lag=12)
autocorTest(residuals(ford.mod11,standardize=T)^2,lag=12)
archTest(residuals(ford.mod11,standardize=T),lag=12)

plot(ford.mod11, ask=F)

# exponential garch models
hp.egarch = garch(hp.s~1, ~egarch(1,1),
leverage=T, trace=F)
hp.egarch
# t-stat on coef
coef(hp.egarch)[5]/sqrt(vcov(hp.egarch)[5,5])

# threshold garch models
hp.tgarch = garch(hp.s~1, ~tgarch(1,1), trace=F)
hp.tgarch
#t-stat on coefficient
coef(hp.tgarch)[5]/sqrt(vcov(hp.tgarch)[5,5])

# power garch
hp.garch = garch(hp.s~1, ~garch(1,1), leverage=T,
trace=F)
hp.garch
# t-stat on leverage effect
coef(hp.garch)[5]/sqrt(vcov(hp.garch)[5,5])

# power garch with d=1
hp.pgarch = garch(hp.s~1, ~pgarch(1,1,1), leverage=T,
trace=F)
hp.pgarch
# t-stat on leverage effect
coef(hp.pgarch)[5]/sqrt(vcov(hp.pgarch)[5,5])

# news impact curves
a0 = hp.tgarch$coef[2]
a1 = hp.tgarch$coef[3]
b1 = hp.tgarch$coef[4]
g1 = hp.tgarch$coef[5]
A = a0 + b1*hp.tgarch$asymp.sd^2
epsilon = seq(-0.21,0.14, length=100)
sigma2.t.TGARCH = A+(a1+g1*(epsilon < 0))*(epsilon^2)

a0 = hp.pgarch$coef[2]
a1 = hp.pgarch$coef[3]
b1 = hp.pgarch$coef[4]
g1 = hp.pgarch$coef[5]
A = (a0 + b1*hp.pgarch$asymp.sd)^2
error = abs(epsilon) + g1*epsilon
sigma2.t.PGARCH = A + 2*sqrt(A)*a1*error + (a1*error)^2

a0 = hp.garch$coef[2]
a1 = hp.garch$coef[3]
b1 = hp.garch$coef[4]
g1 = hp.garch$coef[5]
A = a0 + b1*hp.garch$asymp.sd^2
error = abs(epsilon) + g1*epsilon
sigma2.t.GARCH = A + a1*(error^2)

matplot(cbind(epsilon, epsilon, epsilon),
cbind(sigma2.t.TGARCH,sigma2.t.PGARCH,sigma2.t.GARCH),
type="l")

key(-0.05, 0.0045, lines=list(type="l", lty=1:3),
text=list(c("TGARCH","PGARCH","GARCH")),border=1, adj=1)

# two components model

# persistence parameter and half-life
ford.mod11$coef[3] + ford.mod11$coef[4]
log(0.5)/log(ford.mod11$coef[3] + ford.mod11$coef[4])

# estimate 2 components model pgarch
ford.2comp = garch(ford.s~1, ~pgarch.2comp(2))
summary(ford.2comp)

# 2 components pgarch with leverage effects
garch(ford.s~1, ~two.comp(2,2),leverage=T,trace=F)

# garch in mean

hp.gmean = garch(hp.s~var.in.mean, ~pgarch(1,1,1), leverage=T)
summary(hp.gmean)

# garch with exogenous variables in mean eqn
# capm with garch errors

rvfPlot(ford.s,nyse.s,strip.text="Market Beta",
id.n = 0, hgrid=T, vgrid=T,
xlab="NYSE Returns", ylab="Ford Returns")

ford.beta = garch(ford.s~ma(1)+seriesData(nyse.s),
~garch(1,1), trace=F)
summary(ford.beta)

# garch with exogenous variables in conditional variance

log.abs.ret = log(abs(dell.s-mean(dell.s)))[-1]
# log difference of volume
d.volume = getReturns(dell.v)
# plot abs returns vs volume
rvfPlot(log.abs.ret,d.volume,strip="Scatter Plot",
id.n=0, hgrid=T,vgrid=T,
xlab="% Volume", ylab="Volatility")
# fit garch with volume as exog variable in conditional variance
dell.mod = garch(dell.s~1, ~egarch(1,1)+seriesData(d.volume),
series.start=2)
summary(dell.mod)

# Non Gaussian error distributions

# garch with Student-t errors
ford.mod11.t = garch(ford.s~1, ~garch(1,1), cond.dist="t")
ford.mod11.t
names(ford.mod11.t)
ford.mod11.t$cond.dist
# make qq plot against Student-t: option 13
plot.garch(ford.mod11.t,ask=F,which.plots=12)

# garch with double exponential errors
ford.mod11.dexp = garch(ford.s~1, ~garch(1,1), 
cond.dist="double.exp")
ford.mod11.dexp
names(ford.mod11.dexp)
ford.mod11.dexp$cond.dist
plot.garch(ford.mod11.dexp,ask=F,which.plots=12)
	
# garch with GED errors 
ford.mod11.dexp = garch(ford.s~1, ~garch(1,1),
cond.dist="ged", dist.par=1, dist.est=F)
ford.mod11.dexp
ford.mod11.dexp$cond.dist
plot.garch(ford.mod11.t,ask=F,which.plots=12)

# garch model comparisons

ford.compare = compare.mgarch(ford.mod11,ford.mod11.t)
oldClass(ford.compare)
ford.compare
plot(ford.compare)
plot(ford.compare,qq=T)
	
# garch predictions
hp.pgarch.pred = predict(hp.pgarch,10)
class(hp.pgarch.pred)
names(hp.pgarch.pred)
hp.pgarch.pred