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

intel.daily=read.csv("DailyINTCfrom2000.csv",header=T)
intel.prices=intel.daily$Adj.Close[length(intel.daily$Adj.Close):1]     ## I am rewriting the earlier variables
intel.returns=returns(intel.prices)[-1]

model.intel.garchinM=garchOxFit(formula.mean=~arma(0,0),formula.var=~garch(1,1),series=intel.returns,arch.in.mean=T)


model.intel.egarch=garchOxFit(formula.mean=~arma(0,0),formula.var=~egarch(1,1),series=intel.returns)

#It is a bit more complicated to fit an EGARCH model in R.
#First, an EGARCH(m, s) model in the textbook is an EGACRH(m, s+1) model in R.
#Thus, an EGARCH(1,1) model in R also contains the ARCH parameter.
#Second, the names alpha and beta are interchanged between textbook and R.
#Third, the default option of parameter constraints in R requires the constant term in the volatility equation to be positive.
#This is incorrect because EAGRCH uses logarithm of volatility. Thus, one needs to use “unconstrainted” estimation to fit an EGARCH model.
#This involves editing the “GarchOxModelling.ox” file in the OX/lib directory.
#Finally, THETA1 in R output is the leverage parameter.
#[Alternatively, you may create a separate GarchOxModelling.ox file and garchoxfit-R.txt file to perform the EGARCH estimation in R.]

Strong convergence using numerical derivatives              #We actually fitted an EGARCH(1,2) in the book's notation
Log-likelihood = 584.503
Please wait : Computing the Std Errors ...

 Maximum Likelihood Estimation (Std.Errors based on Second derivatives)
                  Coefficient  Std.Error  t-value  t-prob
Cst(M)              -0.000719  0.0012107  -0.5942  0.5529
Cst(V)               0.000000     593.54     0.00  1.0000
ARCH(Alpha1)         2.226451     2.2623   0.9841  0.3260   # Note that this parameter is not significant so it would seem that EGARCH(1,1) is more appropriate
GARCH(Beta1)         0.991726  0.0033130    299.3  0.0000
EGARCH(Theta1)      -0.006981   0.022285  -0.3133  0.7543   # this is theta in the book's notation
EGARCH(Theta2)       0.188667    0.13711    1.376  0.1700   # this is gamma in the book's notation

No. Observations :       263  No. Parameters  :         6
Mean (Y)         :  -0.00012  Variance (Y)    :   0.00051
Skewness (Y)     :  -0.71426  Kurtosis (Y)    :   7.44730
Log Likelihood   :   584.503

Warning : To avoid numerical problems, the estimated parameter
Cst(V), and its std.Error have been multiplied by 10^4.

source("garchoxfit_R.txt") ## Egarch is also tricky because it blocks the  function. You need to run the source file again to fit other models

## So let us try an EGARCH(1,1):
model.intel.egarch2=garchOxFit(formula.mean=~arma(0,0),formula.var=~egarch(1,0),series=intel.returns)  # Note it is the same output as before - R cannot have s=0

# so the model is:
## log(sigma^2_t) = g(epsilon_{t-1}) /(1-0.992 L)
## g(epsilon_{t-1}) = -0.007 epsilon_{t-1} +0.189 [|epsilon_{t-1}|-sqrt(2/pi)]

###*********************

## TGARCH or GJR model:
source("garchoxfit_R.txt")

model.intel.tgarch=garchOxFit(formula.mean=~arma(0,0),formula.var=~gjr(1,1),series=intel.returns)

Strong convergence using numerical derivatives
Log-likelihood = 632.902
Please wait : Computing the Std Errors ...

 Maximum Likelihood Estimation (Std.Errors based on Second derivatives)
                  Coefficient  Std.Error  t-value  t-prob
Cst(M)              -0.000068  0.0013700 -0.04981  0.9603
Cst(V)               0.324363    0.20996    1.545  0.1236
ARCH(Alpha1)         0.000000   0.070165     0.00  1.0000     #not signiff
GARCH(Beta1)         0.881387   0.078353    11.25  0.0000
GJR(Gamma1)          0.105942   0.063841    1.659  0.0982     # appears marginally nonsigniff

No. Observations :       263  No. Parameters  :         5
Mean (Y)         :  -0.00012  Variance (Y)    :   0.00051
Skewness (Y)     :  -0.71426  Kurtosis (Y)    :   7.44730
Log Likelihood   :   632.902

Warning : To avoid numerical problems, the estimated parameter
Cst(V), and its std.Error have been multiplied by 10^4.

## if we put all param together we get:
## sigma_t^2 =0.0000324363 + (0.000000 + 0.105942 N_1 ) a_{t-1}^2  + 0.881387 sigma_{t-1}^2

### ALTERNATIVE: We can also use the package fGarch from R. The advantage is that using this package we can predict values automatically.

library(fGarch)  # Load the package
n1=garchFit(~garch(1,1),data=intel.returns,trace=F)  # Fit a GARCH(1,1) model
n1   # show the results
predict(n1,5)  # Perform prediction 1 to 5-step ahead.