# Chapter 12 example R code

library(TSA)

## Plot of CREF data:

# Exhibit 12.1
#win.graph(width=4.875, height=2.5,pointsize=8)
data(CREF)
plot(CREF)
t1=477;t2=508
polygon(x=c(t1,t1,t2,t2,t1),
y=c(min(CREF)-10,max(CREF)+10,max(CREF)+10,min(CREF)-10,min(CREF)-10),col='gray')
lines(CREF)

## Plot of returns of CREF data (times 100):

# Exhibit 12.2
r.cref=diff(log(CREF))*100
plot(r.cref)
t1=476;t2=507
polygon(x=c(t1,t1,t2,t2,t1),
y=c(-2,2.7,2.7,-2,-2),col='gray')
lines(r.cref)
abline(h=0)

## ACF and PACF of returns of CREF data:

# Exhibit 12.3
acf(r.cref)

# Exhibit 12.4
pacf(r.cref)

## Nothing out of the ordinary...


## ACF and PACF of absolute returns of CREF data:

# Exhibit 12.5
acf(abs(r.cref))

# Exhibit 12.6
pacf(abs(r.cref))

# We have significant autocorrelations.

## ACF and PACF of squared returns of CREF data:

# Exhibit 12.7
acf(r.cref^2)

# Exhibit 12.8
pacf(r.cref^2)

# We again have significant autocorrelations.

## McLeod-Li test:

# Exhibit 12.9
#win.graph(width=4.875, height=3,pointsize=8)
McLeod.Li.test(y=r.cref)

### Checking for normality of the return series:

# Exhibit 12.10
qqnorm(r.cref)
qqline(r.cref)

shapiro.test(r.cref)


# Implement the Jarque-Bera test for normality in two different ways
# This test is based on the skewness and kurtosis of the sample data:
skewness(r.cref)
kurtosis(r.cref)
length(r.cref)*skewness(r.cref)^2/6
length(r.cref)*kurtosis(r.cref)^2/24
JB=length(r.cref)*(skewness(r.cref)^2/6+kurtosis(r.cref)^2/24)
JB # The Jarque-Bera test statistic
1-pchisq(JB,df=2)

library(tseries)
jarque.bera.test(r.cref)

## All the tests reject the normality of the series.

## Simulated ARCH(1) data:

# Exhibit 12.11
set.seed(1235678)
garch01.sim=garch.sim(alpha=c(.01,.9),n=500)
plot(garch01.sim,type='l',ylab=expression(r[t]),xlab='t')

## Note the apparent volatility clustering.

## Estimating the ARCH(1) model for these simulated data:

m1.sim=garch(x=garch01.sim, order=c(0,1))  # order=c(0,1) specifies an ARCH(1) model
summary(m1.sim)

## Estimates of omega and alpha are: 0.011 and 0.814
## compared to the true values of 0.01 and 0.9 ... not too bad.

## A Shumway and Stoffer example on the GNP data:

library(astsa)

data(gnp)

plot(diff(log(gnp)))  # Plot the return series

acf2(diff(log(gnp))) # How should the model be specified?

gnp.ar1 = sarima(diff(log(gnp)), 1,0,0) # Fit an AR(1) model to the returns

# Look at the plot of the standardized residuals:
# Do we see volatility clustering --> nonconstant variance?

# Let's look at the ACF and PACF of the SQUARED residuals:

acf2(resid(gnp.ar1$fit)^2, 20)

# Maybe some slight autocorrelation (not much, honestly).

# This ACF and PACF of the squared residuals should show an AR(1)-like
# pattern to convince us that the error process is ARCH(1).
# Let's pretend it does...

library(fGarch)  # may need to install 'fGarch' package first

gnp.arch = garchFit(~arma(1,0)+garch(1,0), diff(log(gnp)) ) 

summary(gnp.arch)
# specifies an AR(1) model for the mean process and an ARCH(1) process for the errors

# Estimates are given for the intercept, phi, omega, and alpha.

# Note that the residuals of the model may not be normal.

# Compare these parameter estimates to those from the ordinary AR(1) fit for this series:

print(gnp.ar1)

###################
##
## GARCH examples follow:

## Plot of a simulated GARCH(1,1) series with omega=0.02, alpha=0.05, and beta=0.9:

# Exhibit 12.12
set.seed(1234567)
garch11.sim=garch.sim(alpha=c(0.02,0.05),beta=.9,n=500)
plot(garch11.sim,type='l',ylab=expression(r[t]),xlab='t')

## ACF and PACF of the simulated GARCH(1,1) series:

# Exhibit 12.13
acf(garch11.sim)

# Exhibit 12.14
pacf(garch11.sim)

# Not much autocorrelation ...

## ACF and PACF of the ABSOLUTE simulated GARCH(1,1) series:

# Exhibit 12.15
acf(abs(garch11.sim))

# Exhibit 12.16
pacf(abs(garch11.sim))

# Some substantial autocorrelation ... doesn't look like white noise.

## ACF and PACF of the SQUARED simulated GARCH(1,1) series:

# Exhibit 12.17
acf(garch11.sim^2)

# Exhibit 12.18
pacf(garch11.sim^2)

# Again, substantial autocorrelation ... doesn't look like white noise.

## EACF to try to specify the orders of the GARCH model:

# Exhibit 12.19
eacf((garch11.sim)^2)

# Exhibit 12.20
eacf(abs(garch11.sim))


# Exhibit 12.22
# arima(abs(abs(garch11.sim)),order=c(1,0,1))
# This is fitting an ARMA(1,1) model to the absolute values of the data.
# Just for comparing to the GARCH estimated coefficients found below.

## Fitting the GARCH(2,2) model to the simulated data:

# Exhibit 12.23
g1=garch(garch11.sim,order=c(2,2))
summary(g1)

# Some nonsignificant coefficients!

## Fitting the GARCH(1,1) model to the simulated data:

# Exhibit 12.24
g2=garch(garch11.sim,order=c(1,1))
summary(g2)

AIC(g1)
AIC(g2)

## All in all, the GARCH(1,1) fit seems better.

## The estimated coefficients are omega=0.0075, alpha1= 0.0471 , beta1= 0.935
## compared to the true values of omega=0.02, alpha=0.05, and beta=0.9.  Not too bad.

###########################################
## Fitting a GARCH model to the CREF data:

eacf((r.cref)^2)

# Exhibit 12.21
eacf(abs(r.cref))

# Exhibit 12.25
m1=garch(x=r.cref,order=c(1,1))
summary(m1)

var(r.cref)  # Compare to 0.01633 / (1-0.04414 - 0.91704)

# Exhibit 12.26
plot(residuals(m1),type='h',ylab='standardized residuals')

# Exhibit 12.27
qqnorm(residuals(m1))
qqline(residuals(m1))

# Exhibit 12.28
acf(residuals(m1)^2,na.action=na.omit)

# Exhibit 12.29
gBox(m1,method='squared')

gBox(m1,lags=20, plot=F,method='squared')$pvalue

# Exhibit 12.30
acf(abs(residuals(m1)),na.action=na.omit)
# Exhibit 12.31
gBox(m1,method='absolute')

# Overfitting the GARCH(1,2) model to the CREF returns
m2=garch(x=r.cref,order=c(1,2))
summary(m2,diagnostics=F) 
# The summary is based on the summary.garch function
# in the tseries pacakge. Note the Ljung-Box test from the summary is not
# valid; see the text book. 
AIC(m1)
AIC(m2)

#Further model diagnostic of the fitted GARCH(1,1) model for the CREF returns
shapiro.test(na.omit(residuals(m1)))
jarque.bera.test(na.omit(residuals(m1)))
skewness(na.omit(residuals(m1)))
kurtosis(na.omit(residuals(m1)))

# Exhibit 12.32
plot((fitted(m1)[,1])^2,type='l',ylab='conditional variance',xlab='t')
 

## Dow Jones example:

# library(xts) # not needed?

library(astsa)
data(djia)

djiar=diff(log(djia$Close))[-1] # Calculating the return series
plot(djiar)

# Definite volatility clustering!

acf2(djiar,na.action=na.pass)  # some autocorrelation

acf2(djiar^2,na.action=na.pass)  # lots of autocorrelation

library(fGarch)

# Fitting a AR(1)-GARCH(1,1) model:

djia.g <- garchFit(~arma(1,0)+garch(1,1), data=djiar, cond.dist='std')  # models errors with t distribution...

summary(djia.g)

plot(djia.g)

# Select plots 3 and 2:  note the large conditional variance around the time of the financial crisis of 2008.
# Hit 'Esc' to exit.