#
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
# 
#-----
if( !require('Ecdat') ){
   install.packages('Ecdat', dependencies=TRUE, repos='http://cran.rstudio.com/')
}
library(Ecdat)

if( !require('rugarch') ){
   install.packages('rugarch', dependencies=TRUE, repos='http://cran.rstudio.com/')
}
library(rugarch)

save_plots = FALSE

set.seed(0)

# Exercise 8 EPage 471-472
# 
library(rugarch)
library(Ecdat)
data(SP500,package="Ecdat")
returnBlMon = SP500$r500[1805] 
x = SP500$r500[(1804-2*253+1):1804]
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter14/chap_14_exercise_8_r_black_monday.eps", onefile=FALSE, horizontal=FALSE) }
plot(c(x,returnBlMon))
if( save_plots ){ dev.off() }

# Fit an AR(1)/ARCH(1,1) model:
#
spec = ugarchspec(mean.model=list(armaOrder=c(1, 0)), variance.model=list(garchOrder=c(1, 1)), distribution.model='std')
fit = ugarchfit(data=x, spec=spec)
show(fit)

dfhat = coef(fit)[6]
forecast = ugarchforecast(fit, data=x, n.ahead=1)

meanForecast = as.double(fitted(forecast))
standardDeviation = as.double(sigma(forecast))

# Part (a):
z_score = ( returnBlMon - meanForecast ) / standardDeviation
pt( z_score, dfhat )

# Part (b):
res = forecast@model$modeldata$residuals
res_std = res / forecast@model$modeldata$sigma # standardized residuals
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter14/chap_14_exercise_8_standard_residuals.eps", onefile=FALSE, horizontal=FALSE) }
par(mfrow=c(1,3))
plot(res_std)
acf(res_std)
acf(res_std^2)
par(mfrow=c(1,1))
if( save_plots ){ dev.off() }

# Part (c): Look at fitting a AR(1)/ARCH(1) model:
#
spec = ugarchspec(mean.model=list(armaOrder=c(1, 0)), variance.model=list(garchOrder=c(1, 0)), distribution.model='std')
results_c = ugarchfit(data=x, spec=spec)
res = residuals(results_c)
res_std = res / results_c@fit$sigma # standardized residuals
par(mfrow=c(1,3))
plot(res_std)
acf(res_std)
acf(res_std^2)
par(mfrow=c(1,1))

show(results_c)

# Part (d): Fit a AR(1) model and plot statistics:
#
results_d = arima( x, order=c(1,0,0) )

res = residuals(results_d)
res_std = res / sqrt(results_d$sigma2) # standardized residuals
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter14/chap_14_exercise_8_ar_1_plot.eps", onefile=FALSE, horizontal=FALSE) }
par(mfrow=c(1,3))
plot(res_std)
acf(res_std)
acf(res_std^2)
par(mfrow=c(1,1))
if( save_plots ){ dev.off() }


# E 9 EPage 449
# 
library(rugarch)
library(Ecdat)
data(Irates)
r = as.numeric(log(Irates[,2]))
n = length(r)
lagr = r[1:(n-1)]
diffr = r[2:n] - lagr
spec = ugarchspec(mean.model=list(armaOrder=c(1, 0)), variance.model=list(garchOrder=c(1, 1)), distribution.model='std')
fit = ugarchfit(data=diffr, spec=spec)
show(fit)
plot(fit, which='all')


# E 10 EPage 450:
#
library(rugarch)
library(quantmod)
if( TRUE ){
    getSymbols('^GSPC', from='2005-01-01', to='2014-12-31')
    #write.csv(GSPC, file='../../BookCode/Data/GSPC.csv') # save the data so we can work off line
    sp500 = xts(diff(log(GSPC[, 6]))[-1])
}else{
    GSPC = read.csv('../../BookCode/Data/GSPC.csv')
    GSPC$X = NULL
    sp500 = diff(log(GSPC[, 6]))
}
head(GSPC)

plot(sp500)
y = as.numeric(sp500)

# Part(a):
#
acf(y)
Box.test(y, lag=6, type="Ljung")

# Part(b):
#
acf(y^2)
Box.test(y^2, lag=6, type="Ljung")


# Part(c):
#
y = y[1000:2515] # I could not get ugarchfit to converge for the full range of data 

spec = ugarchspec(mean.model=list(armaOrder=c(1, 0)),  variance.model=list(garchOrder=c(0, 1)), distribution.model='norm')
fit = ugarchfit(data=y, spec=spec) #, solver.control=list(n.restarts=10, maxeval=25000))
show(fit)


# Part (e):
#
spec = ugarchspec(mean.model=list(armaOrder=c(1, 0)), variance.model=list(garchOrder=c(1, 1)), distribution.model='norm')
fit = ugarchfit(data=y, spec=spec)
show(fit)


# Part (f):
#
res = residuals(fit)
res_std = res / fit@fit$sigma # standardized residuals
par(mfrow=c(1,3) )
plot(res_std)
acf(res_std)
acf(res_std^2)
par(mfrow=c(1,1))

Box.test(res_std, lag=20, type='Ljung')
Box.test(res_std^2, lag=20, type='Ljung')


# Part (h):
#
qqnorm(res_std)
qqline(res_std)


# Part (i):
#
plot(fit@fit$sigma)
grid()