[R] Maximum likelihood estimation of ARMA(1,1)-GARCH(1,1)

[Andy Yeh]

Hello

Following some standard textbooks on ARMA(1,1)-GARCH(1,1) (e.g. Ruey
Tsay's Analysis of Financial Time Series), I try to write an R program
to estimate the key parameters of an ARMA(1,1)-GARCH(1,1) model for
Intel's stock returns. For some random reason, I cannot decipher what
is wrong with my R program. The R package fGarch already gives me the
answer, but my customized function does not seem to produce the same
result.

I would like to build an R program that helps estimate the baseline
ARMA(1,1)-GARCH(1,1) model. Then I would like to adapt this baseline
script to fit different GARCH variants (e.g. EGARCH, NGARCH, and
TGARCH). It would be much appreciated if you could provide some
guidance in this case. The code below is the R script for estimating
the 6 parameters of an ARMA(1,1)-GARCH(1,1) model for Intel's stock
returns. At any rate, I would be glad to know your thoughts and
insights. If you have a similar example, please feel free to share
your extant code in R. Many thanks in advance.

Emily



# This R script offers a suite of functions for estimating  the
volatility dynamics based on the standard ARMA(1,1)-GARCH(1,1) model
and its variants.
# The baseline ARMA(1,1) model characterizes the dynamic evolution of
the return generating process.
# The baseline GARCH(1,1) model depicts the the return volatility
dynamics over time.
# We can extend the GARCH(1,1) volatility model to a variety of
alternative specifications to capture the potential asymmetry for a
better comparison:
# GARCH(1,1), EGARCH(1,1),  NGARCH(1,1), and TGARCH(1,1).

options(scipen=10)

intel= read.csv(file="intel.csv")
summary(intel)

raw_data= as.matrix(intel$logret)

library(fGarch)
garchFit(~arma(1,1)+garch(1,1), data=raw_data, trace=FALSE)


negative_log_likelihood_arma11_garch11=
function(theta, data)
{mean =theta[1]
 delta=theta[2]
 gamma=theta[3]
 omega=theta[4]
 alpha=theta[5]
 beta= theta[6]

 r= ts(data)
 n= length(r)

 u= vector(length=n)
 u= ts(u)
 u[1]= r[1]- mean

 for (t in 2:n)
 {u[t]= r[t]- mean- delta*r[t-1]- gamma*u[t-1]}

 h= vector(length=n)
 h= ts(h)
 h[1]= omega/(1-alpha-beta)

 for (t in 2:n)
 {h[t]= omega+ alpha*(u[t-1]^2)+ beta*h[t-1]}

 #return(-sum(dnorm(u[2:n], mean=mean, sd=sqrt(h[2:n]), log=TRUE)))
 pi=3.141592653589793238462643383279502884197169399375105820974944592
 return(-sum(-0.5*log(2*pi) -0.5*log(h[2:n]) -0.5*(u[2:n]^2)/h[2:n]))
}


#theta0=c(0, +0.78, -0.79, +0.0000018, +0.06, +0.93, 0.01)
theta0=rep(0.01,6)
negative_log_likelihood_arma11_garch11(theta=theta0, data=raw_data)


alpha= proc.time()
maximum_likelihood_fit_arma11_garch11=
nlm(negative_log_likelihood_arma11_garch11,
    p=theta0,
    data=raw_data,
    hessian=TRUE,
    iterlim=500)
#optim(theta0,
#      negative_log_likelihood_arma11_garch11,
#      data=raw_data,
#      method="L-BFGS-B",
#      upper=c(+0.999999999999,+0.999999999999,+0.999999999999,0.999999999999,0.999999999999,0.999999999999),
#      lower=c(-0.999999999999,-0.999999999999,-0.999999999999,0.000000000001,0.000000000001,0.000000000001),
#      hessian=TRUE)

# We record the end time and calculate the total runtime for the above work.
omega= proc.time()
runtime= omega-alpha
zhours = floor(runtime/60/60)
zminutes=floor(runtime/60- zhours*60)
zseconds=floor(runtime- zhours*60*60- zminutes*60)
print(paste("It takes ",zhours,"hour(s)", zminutes," minute(s) ","and
", zseconds,"second(s) to finish running this R program",sep=""))


maximum_likelihood_fit_arma11_garch11

sqrt(diag(solve(maximum_likelihood_fit_arma11_garch11$hessian)))