/* This code estimates the NIG-ACD model in Wilhelmsson (2007) "Value-at-risk with time-varying variance, skewness and kurtosis, The NIG-ACD model" The Econometrics Journal 2009 
This version 3 December 2007
*/

// Note parameter restrictions are not enforced by this code.

load data[] =d:\sp500all.txt; // Change this to your dataset
z=data;

library cml;
cmlset;
_cml_gradorder=3;
// Polynomial approximation (from Abromiwitz and Stegun) of the modified Besselfunction of order three with index 1. This function does not exist in Gauss.
proc (1)=besselk(z);
local y,ans,nobs,i,x;
nobs=rows(z);
ans=zeros(nobs,1);
for i(1,nobs,1);
x=z[i];
if x<=2;
  y = x*x/4.0;
ans[i]=ln(x/2.0)*mbesseli1(x)+(1.0/x)*(1.0+y*(0.15443144+y*(-0.67278579+y*(-0.18156897+y*(-0.1919402e-1+y*(-0.110404e-2+y*(-0.4686e-4)))))));
else; 
 y = 2.0/x;
 ans[i] = (exp(-x)/sqrt(x))*(1.25331414+y*(0.23498619+y*(-0.3655620e-1+y*(0.1504268e-1+y*(-0.780353e-2 +y*(0.325614e-2+y*(-0.68245e-3)))))));
endif;
endfor;
retp(ans);
endp;

// The log likelihood function of the NIG-ACD model 
  

proc BNgarch(theta,data);

local LL,u,rhotilde,rho,vgamma,beta,alpha,hf,tal,delett,qform,lltemp,nobs,a,b,c,v,tau,my,d,e,f,g,gg,y;
a=theta[2]; 
b=theta[3]; 
c=theta[4]; 
v=theta[6]; 
tau=theta[5]; 
my=theta[1]; 
d=theta[10]; 
e=theta[11]; 
f=theta[7]; 
g=theta[8]; 
gg=theta[9];
y=data;
nobs=rows(y)-1;
vgamma=ones(nobs,1); 
beta=zeros(nobs,1); 
hf=ones(nobs,1); 
u=y[1:nobs,1]; 
rho=zeros(nobs,1);
LL=zeros(nobs+1,1);

for i(1,nobs-1,1);
 
  rho[i+1,1] = -(exp(d + e*u[i,1])-1)/(-exp(d + e*u[i,1])-1);
  vgamma[i+1] = exp(f + g*u[i]^2+gg*ln(vgamma[i]));
  hf[i+1] = (c + a*(abs(u[i])-tau*u[i])^v + b*hf[i]^v)^(1/v);
  u[i+1] = y[i+1,1] -my - hf[i+1]*sqrt(vgamma[i+1])*rho[i+1]; 

endfor;
beta=rho.*vgamma./(sqrt(1-rho.^2));

retp(0|(2*ln(sqrt(beta.^2+vgamma.^2))-ln(hf)-1.5*ln(vgamma)+vgamma+beta.*(sqrt(beta.^2+vgamma.^2).*u./
(hf.*vgamma.^1.5)+beta./vgamma)-ln(((1+(sqrt(beta.^2+vgamma.^2).*u./(hf.*vgamma.^1.5)+beta./vgamma).^2).^0.5))+
ln(besselk(sqrt(beta.^2+vgamma.^2).*((1+(sqrt(beta.^2+vgamma.^2).*u./(hf.*vgamma.^1.5)+beta./vgamma).^2).^0.5)))-ln(pi))); 

endp;

__weight = (rows(z)/(rows(z)-1))*ones(rows(z),1);
__weight[1,1] = 0;   
theta0={0.0842,0.07,0.93,0.009,0.89,1.0,0.87,-0.18,0.70,-0.06,0.21};  // startingvalues for the parameters
       
{b,f,g,vcv,ret } = cml( z, 0, &bngarch, theta0 );
b;