' Estimates ARMA(p,q)-GARCH(1,1)-Gaussian
' Define the AR and MA lag lengths by adjusting R and M 
' Workfile should be opened and prg file must be in default path  

!R = 3 ' AR lag length
!M=3 ' MA lag length

' Store AR and MA lags
%arma = ""

for !i=1 to !R
   %arma = %arma +" ar(" +@str(!i) +")"
next

for !j=1 to !M
%arma = %arma +" ma(" +@str(!j) +")"
next


'Log return series
series y=dlog(ise)

' declare coef vectors to use in MLE
' mu, phi and theta are mean eqn parameters (constant, AR lags, MA lags)
' omega, alpha, beta are GARCH parameters (constant, resid, sigma)
coef(1) mu = 0.005
coef(!R) phi =0.005 ' AR parameters
coef(!M) theta = 0.005 ' MA parameters
coef(1) e_omega = log(0.005) ' omega=exp(e_omega) (to satisfy omega>0)
coef(1) e_alpha = log(0.005) ' alpha=exp(e_alpha) (to satisfy alpha>0)
coef(1) e_beta = log(0.005) ' beta=exp(e_beta) (to satisfy beta>0)

' Initial guess for ARMA coeffs from OLS
equation eq_ARMA.ls y c {%arma}
mu(1)=eq_ARMA.c(1)

for !i=1 to !R
phi(!i)=eq_ARMA.c(1+!i)
next

for !i=1 to !M
theta(!i)=eq_ARMA.c(1+!R+!i)
next

' Initial guess for variance coeffs (adhoc)
e_alpha(1)=log(0.15)
e_beta(1)=log(0.60)
e_omega(1)=log(0.75*eq_ARMA.@se^2) ' To make theoretical uncond variance to sample error variance

' Use Eviews starting values
mu(1)=0.00072
phi=0.005
theta=0.005
e_alpha(1)=log(0.15)
e_beta(1)=log(0.60)
e_omega(1)=log(0.00028) 

' Set presample values for residuals and variance
smpl 1 !R  
series sig2 = eq_ARMA.@se^2
series res=0
series resma = 0
smpl @all

' Set up log-likelihood function
logl GARCH_LLF
GARCH_LLF.append @logl LLF
GARCH_LLF.append res=y-mu(1)

GARCH_LLF.append resma = res ' To make resma = res - phi(-1)*res(-1) - ... - theta(-1)*resma(-1)
for !i=1 to !R
GARCH_LLF.append resma = resma - phi(!i)*res(-!i)
next
for !j=1 to !M
GARCH_LLF.append resma = resma - theta(!j)*resma(-!j)
next

GARCH_LLF.append sig2= exp(e_omega(1)) + exp(e_alpha(1))*resma(-1)^2 +exp(e_beta(1))*sig2(-1)
GARCH_LLF.append z=resma/@sqrt(sig2)
GARCH_LLF.append LLF=log(@dnorm(z))-log(sig2)/2


' Estimate and display results
smpl !R+1 @last  
GARCH_LLF.ml(showopts, b, m=5000, c=1e-5) 'b means BHHH (default marquardt), m max iter, c convergence criteria
show GARCH_LLF.output

' Estimate the same model by EViews 
equation eq_EViews
eq_EViews.arch(1,1,showopts, z, b, m=5000, c=1e-5) y c {%arma} 
show eq_EViews