```# -------------------------------------------------------------------
# Lecture 4: Modelling volatility of S&P returns

# Required libraries: zoo, lmtest, TSA, rgarch, moments
rm(list=ls())
source("http://klein.uk/R/myfunctions.R")
ls()
# -------------------------------------------------------------------

## S&p500 1988-2010
# When we study stock returns, the main target is to model volatility as stock price may follow a random walk
str(sp)

date <- as.Date(sp\$Date, format="%d/%m/%Y")
library(zoo)

# As there is an extreme value of volatility in 1987, we select the sample from 1988 to 2010 to
# estimate the Garch model.

## Step 1:
# Plot dlsp
plot(dlsp)
# Do unit root test
dlsp.v <- as.vector(dlsp)
# -> DLSP is stationary

## Step2:
# Fit a mean model and check whether the residuals have arch effects.
# First look at correlogram of DLSP
par(mfrow=c(2,1))
acf(na.omit(dlsp.v)); pacf(na.omit(dlsp.v))
# Fit arma model: dlsp = c ar(1) ar(2) ma(1) ma(2)
arima202 <- arima(dlsp, order=c(2,0,2))
library(lmtest)
coeftest(arima202)
# Correlogram of standardized residuals
nNA <- which(is.na(arima202\$resid)==F)
# install.packages("TSA", dependencies=T)
library(TSA)
acf(rstandard(arima202)[nNA]); pacf(rstandard(arima202)[nNA])
# -> Still correlations
# Correlogram of sqared standardized residuals
acf(rstandard(arima202)[nNA]^2); pacf(rstandard(arima202)[nNA]^2)
# -> Still arch effects

## Step 3:
# Fit one arch model to clean arch effects, start with the longest lags
# dlsp c ar(1) ar(2) ma(1) ma(2) Arch(9)

# manually install rgarch and dependencies as described on
# http://klein.uk/Econometrics2.html
library(rgarch)

# ?ugarchspec
spec <- ugarchspec(variance.model = list(model = "fGARCH", submodel="GARCH",
garchOrder = c(9,0)),
mean.model = list(armaOrder = c(2,2), include.mean=T))
# ?ugarchfit
sgarch.fit <- ugarchfit(data=dlsp.v, spec = spec)
sgarch.fit

# Correlogram of residuals
nNA <- which(is.na(sgarch.fit@fit\$resid)==F)
stdresid <- sgarch.fit@fit\$residuals/sgarch.fit@fit\$sigma
acf(stdresid[nNA]); pacf(stdresid[nNA])
ljung.box.test.1(stdresid[nNA], seq(1,15,1))
# -> Still correlations

# Correlogram of sqared residuals
acf(stdresid[nNA]^2); pacf(stdresid[nNA]^2)
ljung.box.test.1(stdresid[nNA]^2, seq(1,15,1))
# -> No arch effects

## Step 4:
# So we try different models to clean correlations among residuals and arch effects
# and find

# dlsp = ar(1) ma(1) ma(12) Garch(2,1)

spec <- ugarchspec(variance.model = list(model = "fGARCH", submodel="GARCH",
garchOrder = c(2,1)),
mean.model = list(armaOrder = c(1,12), include.mean=T),
fixed.pars = list(ma2=0, ma3=0, ma4=0, ma5=0, ma6=0, ma7=0, ma8=0, ma9=0, ma10=0, ma11=0))
sgarch.fit <- ugarchfit(data=dlsp, spec = spec)
sgarch.fit

# Has no problem with correlation in standardized residuals.
# And all coefficients are significant.
nNA <- which(is.na(sgarch.fit@fit\$resid)==F)
stdresid <- sgarch.fit@fit\$residuals/sgarch.fit@fit\$sigma
acf(stdresid[nNA]); pacf(stdresid[nNA])
ljung.box.test.1(stdresid[nNA], seq(1,15,1))

## Step 5:
# Check normality of standardized residuals and do forecasting
# Huge improvements in distributions in terms of Jarque-Bera statistics
library(moments)
jarque.test(stdresid)
# This is the original distribution
jarque.test(dlsp.v)

# Static forecast: see the conditional variance changes
#1
?ugarchforecast

## Step 6:
# However there may be still arch mean effects and leverage effects.
# So to improve our model, we may try TARCH model, IGARCH model , Garch_mean model,
# EGARCH_mean model.

# 1) TARCH model [OPTIONAL]
?ugarchspec
spec <- ugarchspec(variance.model = list(model = "fGARCH", submodel="TGARCH",
garchOrder = c(2,1)),
mean.model = list(armaOrder = c(1,12), include.mean=T),
fixed.pars = list(ma2=0, ma3=0, ma4=0, ma5=0, ma6=0, ma7=0, ma8=0, ma9=0, ma10=0, ma11=0))
sgarch.fit <- ugarchfit(data=dlsp.v, spec = spec)
sgarch.fit

# 2) IGARCH model [OPTIONAL]
spec <- ugarchspec(variance.model = list(model = "iGARCH", garchOrder = c(2,1)),
mean.model = list(armaOrder = c(1,12), include.mean=T),
fixed.pars = list(ma2=0, ma3=0, ma4=0, ma5=0, ma6=0, ma7=0, ma8=0, ma9=0, ma10=0, ma11=0))
sgarch.fit <- ugarchfit(data=dlsp.v, spec = spec)
sgarch.fit
# has problems with correlation in residuals

# 3) EGRACH model [OPTIONAL]
# After many trials to capture the arch mean effects and leverage effects,
# We obtain a good model: dlsp = ar(1) ma(1) ma(12) EGARCH(2,1,1)
spec <- ugarchspec(variance.model = list(model = "eGARCH", garchOrder = c(2,1)),
mean.model = list(armaOrder = c(1,12), include.mean=T),
fixed.pars = list(ma2=0, ma3=0, ma4=0, ma5=0, ma6=0, ma7=0, ma8=0, ma9=0, ma10=0, ma11=0))
sgarch.fit <- ugarchfit(data=dlsp.v, spec = spec)
sgarch.fit

## Interpretation of coefs!
# The negative news RESID(-1) will increase the conditional volatility (ln(Garch) by 0.092359-0.000866.(-
# c(8)+c(6)). The positive news RESID(-1) will decrease the conditional volatility (ln(Garch) by -0.092359-
# 0.000866.( c(8)+c(6))

## 4) Garch in mean [USEFUL]
spec <- ugarchspec(variance.model = list(model = "fGARCH", submodel="GARCH",
garchOrder = c(2,1)),
mean.model = list(armaOrder = c(1,12), include.mean=T, garchInMean=T),
fixed.pars = list(ma2=0, ma3=0, ma4=0, ma5=0, ma6=0, ma7=0, ma8=0, ma9=0, ma10=0, ma11=0))
sgarch.fit <- ugarchfit(data=dlsp, spec = spec)
sgarch.fit
# Interpretation: There is arch mean effects showing during the volatile periods, the risk premium rises as
# the coefficient of garch is 0.066680 in the mean model. This is a time varying risk premium.

# No correlation in standardized residuals:
nNA <- which(is.na(sgarch.fit@fit\$resid)==F)
stdresid <- sgarch.fit@fit\$residuals/sgarch.fit@fit\$sigma
ljung.box.test.1(stdresid[nNA], seq(1,15,1))
# No arch effects in standardized residuals:
ljung.box.test.1(stdresid[nNA]^2, seq(1,15,1))

# Improvements in distributions in terms of Jarque-Bera statistics
jarque.test(stdresid)

# Static forecast: see the conditional variance changes
?ugarchforecast