library(fSeries)
source("garchoxfit_R.txt")

monthly.sp500=read.csv("Monthlysp500.csv",header=T)
monthly.sp50.prices= monthly.sp500$Adj.Close[length(monthly.sp500$Adj.Close):1]
monthlysp500.returns=returns(monthly.sp50.prices)[-1]

acf(monthlysp500.returns)
pacf(monthlysp500.returns)

## I am going to use AR(5)  for the return
## Why garch(1,1) ? order for garch is hard to identify. NO particular reason except that the book uses the same order.

sp500.model1=garchOxFit(formula.mean=~arma(5,0),formula.var=~garch(1,1),series=monthlysp500.returns)
#                  Coefficient  Std.Error  t-value  t-prob
#Cst(M)               0.006157  0.0015488    3.975  0.0001
#AR(1)                0.018841   0.040518   0.4650  0.6421
#AR(2)               -0.046283   0.040166   -1.152  0.2496
#AR(3)                0.010415   0.039289   0.2651  0.7910
#AR(4)                0.025837   0.040014   0.6457  0.5187
#AR(5)                0.081760   0.039317    2.080  0.0379
#Cst(V)               0.771007    0.31817    2.423  0.0156
#
#ARCH(Alpha1)         0.092019   0.024360    3.777  0.0002
#GARCH(Beta1)         0.866546   0.027808    31.16  0.0000

##
sp500.model2=garchOxFit(formula.mean=~arma(5,0),formula.var=~garch(1,1),series=monthlysp500.returns,cond.dist = "t")
#                  Coefficient  Std.Error  t-value  t-prob
#Cst(M)               0.007754  0.0015226    5.093  0.0000
#AR(1)               -0.005713   0.038764  -0.1474  0.8829
#AR(2)               -0.039462   0.038882   -1.015  0.3105
#AR(3)                0.009739   0.038608   0.2523  0.8009
#AR(4)                0.032466   0.038070   0.8528  0.3941
#AR(5)                0.101868   0.037920    2.686  0.0074
#Cst(V)               1.224162    0.54803    2.234  0.0258
#ARCH(Alpha1)         0.099817   0.031166    3.203  0.0014
#GARCH(Beta1)         0.830252   0.046042    18.03  0.0000
#Student(DF)          7.430838     1.9605    3.790  0.0002

m1=arima(monthlysp500.returns,order=c(5,0,0))
# m1
#
#Call:
#arima(x = monthlysp500.returns, order = c(5, 0, 0))
#
#Coefficients:
#         ar1      ar2     ar3     ar4     ar5  intercept
#      0.0258  -0.0486  0.0199  0.0087  0.1014     0.0061
#s.e.  0.0375   0.0377  0.0377  0.0377  0.0377     0.0017
#
#sigma^2 estimated as 0.001648:  log likelihood = 1253.22,  aic = -2492.44


pacf(m1$resid^2)
acf(m1$resid^2)
m1.resid=arima(m1$resid,order=c(1,0,1))        ## quite different

## fit an IGARCH
sp500.model3= garchOxFit(formula.mean=~arma(5,0),formula.var=~igarch(1,1),series=monthlysp500.returns)

## GARCH-M
sp500.model4=garchOxFit(formula.mean=~arma(5,0),formula.var=~garch(1,1),series=monthlysp500.returns,arch.in.mean=T)