I want to estimate a GARCH(1,1) proces manually using maximum likelihood function which is also related to my previous post.
I wrote the below code:
library(fGarch)
#generate returns
set.seed(1)
garch11<-garchSpec(model = list())
y<-garchSim(garch11, n = 1000)
n<-length(y)
sigma2<- vector(length = n)
#estimating the parameters manually by mle
LL <- function(theta)
{
sigma2[1]<-theta[1]/(1-theta[2]-theta[3])
#sigma2[1]<-var(y)
for(i in 1:n)
{
sigma2[i+1]<-theta[1]+theta[2]*y[i]^2+theta[3]*sigma2[i]
}
R =dnorm(y[1:n], 0, sqrt(sigma2[1:n]), log = TRUE)
L=-sum(R)
return(L)
}
fit1<-optim(theta<-c(0.1,0.,0.1),LL)
fit2<-optim(theta<-c(0.0001,0.5,0.1),LL)
> fit1
$par
[1] 0.000007919175 0.139449431273 0.190760467531
$value
[1] -4264.499
$counts
function gradient
250 NA
$convergence
[1] 0
$message
NULL
> fit2
$par
[1] 0.000001281119 0.100213870720 0.790732098375
$value
[1] -4275.474
$counts
function gradient
314 NA
$convergence
[1] 0
$message
NULL
I also wanted to do a manual MLE once and found the mle()
function very unstable, due to the sensitivity to the starting value choice. Trying the solnp()
function from the Rsolnp
package worked for me - here's the code:
#data <- "your data set"
library(Rsolnp)
solfun<-function(vec)
{
alpha<-vec[1];
beta<-vec[2];
omega<-vec[3];
mu<-vec[4]; #only if you want to include a mean in the GARCH model, else set mu=0
sig0<-vec[5]; #starting value estimation can be replaced by sd(data)
LL_ret<-1:(length(data));
sigs<-1:(length(data)+1);
sigs[1]<-sig0^2 ;
for(i in 2:(length(data)+1)) (sigs[i] <- (omega+alpha*(data[i-1]-mu)^2+beta*sigs[i-1])) ;
for(i in 1:(length(data))) (LL_ret[i] <- (-0.5*log(2*pi*sigs[i])-((data[i]-mu)^2)/(2*sigs[i]))) ;
return(sum(-LL_ret)) ;
} #negative log likelihood that we want to minimize
#optimization:
sol<-solnp(pars=c(alpha=0.1, beta=0.9, omega=0, mu=mean(data),sig0=sd(data)),
fun=solfun,LB=c(0,0.7,0,-1,0),UB=c(0.3,1,1,1,1)) #mu and sig0 only if you use them in solfun
#plugging the estimated parameters into function for sigma^2 estimates:
LL2<-function(alpha, beta, omega, mu, sig0)
{
LL_ret<-1:(length(data));
sigs<-1:(length(data)+1);
sigs[1]<-sig0^2 ;
for(i in 2:(length(data)+1)) (sigs[i] <- (omega+alpha*(data[i-1]-mu)^2+beta*sigs[i-1])) ;
for(i in 1:(length(data))) (LL_ret[i] <- (-0.5*log(2*pi*sigs[i])-((data[i]-mu)^2)/(2*sigs[i]))) ;
return(list(LL_ret,sigs)) ;
}
L2<-LL2(alpha=sol$pars[[1]],beta=sol$pars[[2]],
omega=sol$pars[[3]],mu=sol$pars[[4]],sig0=sol$pars[[5]])
#Summarizing results of log likelihood, parameters and sigma estimates:
GG1<-list(sigma=sqrt(L2[[2]]),llh=L2[[1]],llhs=sum(L2[[1]]),coef=sol$pars)
If you don't want to estimate the mean and starting value, you just have to remove those parts. Also, you'll obviously have to adjust the parameter bounds in the optimization, to fit your data (in my case, it was financial return data). It doesn't have to be as exact, as for the mle()
function, though.