/*
** GPE Application Module: GARCH1.GPE
** nonlinear maximum likelihood estimation
** for GARCH(1,1), GARCH-M(1,1) based on normal and
** non-normal distributions (e.g., t and GED distributions). 
**
** (C) Copyright 2001 by Applied Data Associates
** All Rights Reserved.
**
** This Software Product is PROPRIETARY SOURCE CODE OF APPLIED
** DATA ASSOCIATES. This File Header must accompany all files
** using any portion, in whole or in part, of this Source Code.
** This Software Product is designed to be used with GPE2 and
** GAUSS. If you use this Source Code for research and development,
** a proper reference is required. If you wish to distribute any 
** portion of the proprietary Source Code, in whole or in part, you 
** must first obtain written permission from Applied Data Associates.
**
** Consider the linear regression equation y = x b + e (GARCH) or
** y = x b + c h + e (GARCH-M), and the conditional variance equation 
** follows Bollerslev's first-order form: h = a0 + a1 h(-1) + d1 e(-1)^2.
** The predefined procedure ngarch1, tgarch1, and egarch1 are the log-
** likelihood function for GARCH(1,1) model based on normal, t-, and 
** GED distributions, respectively. Similarly, ngarchm1, tgarchm1, and 
** egarchm1 are the log-likelihood function for GARCH-M(1,1) model based 
** on normal, t-, and GED distributions, respectively.
** To estimate a GARCH(1,1) model, just call estimate procedure with 
** ngarch1 log-likelihood function as follows:
** 
** ==> call estimate(&ngarch1,z);
**
** Similarly, to estimate GARCH-M(1,1) model,
**
** ==> call estimate(&ngarchm1,z);
**
** The data matrix z =[y,x] is arranged so that the left-hand-side 
** dependent variable of the linear regression equation is in the 
** first column, followed by the right-hand-side explanatory variables 
** including constant. Based on non-normal distribution, tgarch1 
** (assuming t-distribution) and egarch1 (assuming GED distribution) 
** may be used instead. 
**
** All these predefined log-likelihood functions are vector-valued 
** component functions, therefore the control variable _nlopt should 
** be set to 2 before calling estimate procedure as follows:
**
** _nlopt = 2;
**
** The predefined first-order GARCH and GARCH-M models are for the linear 
** regression equation only. It is straightforward to extend the model 
** for nonlinear regression equation and for higher order of GARCH or 
** GARCH-M strucutre. The application module GARCH1.GPE is designed that 
** it can be easily modified and extended.
*/ 

/*
** Simple GARCH(1,1)
** Y=Xb+e, rows(b)=k
** h=b[k+3]+b[k+2]*e(-1)^2+b[k+1]*h(-1)
*/
proc (2) = garch1(x,b);
    local n,k,i,e,h,me2;
    n=rows(x);     @ number of observations @
    k=cols(x)-1;   @ number of RHS regressors @
    e=x[.,1]-x[.,2:k+1]*b[1:k];
    
    /* data initialization */
    @ me2=meanc(e^2); @           @ -with sample mean @
    me2=b[k+3]/(1-b[k+2]-b[k+1]); @ -with expected value @
/* 
    @ recursive h: slower method using loop @
    h=ones(n,1);
    h[1]=b[k+3]+b[k+2]*me2+b[k+1]*me2;
    i=2; 
    do until i>n;
        h[i]=b[k+3]+b[k+2]*(e[i-1]^2)+b[k+1]*h[i-1];
        i=i+1;
    endo;
*/
    h=b[k+3]+b[k+2]*missrv(lagn(e^2,1),me2);
    h=recserar(me2|h,me2,b[k+1]);
    h=h[2:n+1];
    
    retp(e,h);
endp;

/*
** Simple GARCH-M(1,1)
** Y=Xb+b[k+1]*h+e, rows(b)=k
** h=b[k+4]+b[k+3]*e(-1)^2+b[k+2]*h(-1)
*/
proc (2) = garchm1(x,b);
    local n,k,i,e,h,me2;
    n=rows(x);    @ number of observations @
    k=cols(x)-1;  @ number of regressors @
    e=zeros(n,1);
    h=ones(n,1);

    @ data initialization with E(e^2)=E(h) @
    me2=b[k+4]/(1-b[k+3]-b[k+2]);
    e[1]=x[1,1]-x[1,2:k+1]*b[1:k]-b[k+1]*me2;
    h[1]=b[k+4]+b[k+3]*me2+b[k+2]*me2;
    
    i=2; @ e depends on recursive h @
    do until i>n;
        h[i]=b[k+4]+b[k+3]*(e[i-1]^2)+b[k+2]*h[i-1];
        e[i]=x[i,1]-x[i,2:k+1]*b[1:k]-b[k+1]*h[i];
        i=i+1;
    endo;
    retp(e,h);
endp;

/*
** NGARCH(1,1) log-likelihood function
** GARCH(1,1) based on normal density
*/
proc ngarch1(x,b);
    local e,h,ll;
    {e,h}=garch1(x,b);    
    h=(h.>0).*h+(h.<=0).*(1.0e-8); @ make sure h.>0 @
    /*
    if h>0;
        ll=-0.5*(ln(2*pi)+ln(h)+(e^2)./h);
    else; 
        ll=-1e+15; 
    endif;
    */
    retp(-0.5*(ln(2*pi)+ln(h)+(e^2)./h));
endp;

/*
** NGARCH-M(1,1) log-likelihood function
** GARCH-M(1,1) based on normal density
*/
proc ngarchm1(x,b);
    local e,h,ll;
    {e,h}=garchm1(x,b);
    h=(h.>0).*h+(h.<=0).*(1.0e-8); @ make sure h.>0 @
    /*
    if h>0;
        ll=-0.5*(ln(2*pi)+ln(h)+(e^2)./h);
    else; 
        ll=-1e+15; 
    endif;
    */
    retp(-0.5*(ln(2*pi)+ln(h)+(e^2)./h));
endp;

/*
** TGARCH(1,1) log-likelihood function
** based on t distribution
** b[k]=d.f. parameter of t-distribution
** b[k] is the last estimated parameter
*/
proc tgarch1(x,b);
    local e,h,k,df,ll;
    k=rows(b);
    df=b[k];
    {e,h}=garch1(x,b[1:k-1]);
    h=(h.>0).*h+(h.<=0).*(1.0e-8); @ make sure h.>0 @
    ll=lnfact((df-1)/2)-lnfact((df-2)/2)-0.5*ln(df-2)-0.5*ln(pi);
    /*
    if h>0;
        ll=ll-0.5*(ln(h)+(df+1)*ln(1+(e^2)./(h.*(df-2))));
    else; 
        ll=ll-1e+15; 
    endif;
    */
    retp(ll-0.5*(ln(h)+(df+1)*ln(1+(e^2)./(h.*(df-2)))));
endp;

/*
** TGARCH-M(1,1) log-likelihood function
** based on t distribution
** b[k]=d.f. parameter of t-distribution
** b[k] is the last estimated parameter
*/
proc tgarchm1(x,b);
    local e,h,k,df,ll;
    k=rows(b);
    df=b[k];
    {e,h}=garchm1(x,b[1:k-1]);
    h=(h.>0).*h+(h.<=0).*(1.0e-8); @ make sure h.>0 @
    ll=lnfact((df-1)/2)-lnfact((df-2)/2)-0.5*ln(df-2)-0.5*ln(pi);
    /*
    if h>0;
        ll=ll-0.5*(ln(h)+(df+1)*ln(1+(e^2)./(h.*(df-2))));
    else; 
        ll=ll-1e+15; 
    endif;
    */
    retp(ll-0.5*(ln(h)+(df+1)*ln(1+(e^2)./(h.*(df-2)))));
endp;

/*
** EGARCH(1,1) log-likelihood function
** based on generalized error distribution
** b[k]>0: thickness of the tails of the GED
** b[k]=2: a standard normal density
** b[k]<2: thicker tails than the normal density
** b[k]>2: thinner tails than the normal density
** b[k] is the last estimated parameter
*/
proc egarch1(x,b);
    local e,h,k,v,l,ll;

    {e,h}=garch1(x,b);
    k=rows(b);
    v=b[k];
    l=sqrt((2^(-2/v))*(gamma(1/v)./gamma(3/v)));      
    ll=ln(v/l)-(1/v+1)*ln(2)-ln(gamma(1/v))
       -0.5*(abs(e./(l*sqrt(h)))^v+ln(h));       
    retp(ll);
endp;

/*
** EGARCH-M(1,1) log-likelihood function
** based on generalized error distribution
** b[k]>0: thickness of the tails of the GED
** b[k]=2: a standard normal density
** b[k]<2: thicker tails than the normal density
** b[k]>2: thinner tails than the normal density
** b[k] is the last estimated parameter
*/
proc egarchm1(x,b);
    local e,h,k,v,l,ll;

    {e,h}=garchm1(x,b);
    k=rows(b);
    v=b[k];
    l=sqrt((2^(-2/v))*(gamma(1/v)./gamma(3/v)));      
    ll=ln(v/l)-(1/v+1)*ln(2)-ln(gamma(1/v))
       -0.5*(abs(e./(l*sqrt(h)))^v+ln(h));       
    retp(ll);
endp;