/*##############################################################################
#                                                                              #
#                 VaR by Breno de Andrade Pinheiro Néri                        #
#                                                                              #
#            brenoneri@gmail.com       www.fgv.br/aluno/bneri                  #
#                                                                              #
#                                                    August, 2005 - Version 7  #
#                                                                              #
################################################################################
#                                                                              #
# This code estimates 1-day-ahead-forecast VaR time series using three         #
# methodologies, namely:                                                       #
# 1- RiskMetrics model developed by J.P. Morgan, 1996;                         #
# 2- GARCH(1, 1) with Gaussian errors;                                         #
# 3- APARCH(1, 1) with Skewed Student-t errors.                                #
#                                                                              #
# You enter a time series vector (vY) in a file called vY.mat, with the        #
# observations in one single column, and the code saves a matrix in a file     #
# called mVaR.mat, with the 3 VaR methodologies in the columns and the 1-day   #
# VaR forecast in the rows.                                                    #
#                                                                              #
# Note that each forecasts are estimated based on the last 'cWin' observations.#
# Hence, if rows(vY) < cWin, you will have an error. If you do not have enough #
# observations, set up a lower temporal window size (cWin).                    #
# Note also that rows(mVaR) = rows(vY) - cWin + 1.                             #
#                                                                              #
# Due to the two constrained nonlinear optimization for each observation, this #
# code is extremely computational intensive!                                   #
#                                                                              #
# Defaults: tau=.01 and cWin=250, due to regulatory obligation (1996 Amendment #
# to Basle Capital Accord).                                                    #
#                                                                              #
################################################################################
#                                                                              #
# You may use this code only if you accept the conditions:                     #
# (1) I am not liable for any problem caused by this code (i.e. you use it at  #
# your own risk);                                                              #
# (2) You must give me credit in your papers where this code has been used.    #
#                                                                              #
##############################################################################*/

#import<packages/Garch40/garch>

static decl s_vE, s_obj;

// GARCH11
garch(const vP)
{
    decl vSigma2=s_vE[0]^2, i;

    for(i=0; i<rows(s_vE); ++i) vSigma2 |= vP[0] + vP[1]*s_vE[i]^2 + vP[2]*vSigma2[i];
    return(vSigma2);
}

// APARCH11
aparch(const vP)
{
    decl sigmad=fabs(s_vE[0])^vP[4], i;

    for(i=0; i<rows(s_vE); ++i) sigmad |= vP[0] + vP[1]*(fabs(s_vE[i])-vP[2]*s_vE[i])^vP[4] + vP[3]*sigmad[i];
    return(sigmad.^(2/vP[4]));
}

// Joint Likelihoods
LG(const vP, const dP, const vScore, const mHessian) {return(dP[0]=1e-2*sumc(s_obj.GaussLik(s_vE, garch(vP)[:rows(s_vE)-1])));}
LST(const vP, const dP, const vScore, const mHessian) {return(dP[0]=1e-5*sumc(s_obj.SkStudentLik(s_vE, aparch(vP[0:4])[:rows(s_vE)-1], log(vP[6]), vP[5])));}

main()
{
    decl vY=loadmat("vY.mat", 1), cT=rows(vY), tau=.01, cWin=250, mVaR=zeros(cT-cWin+1,3), quantnorm=quann(tau), i, j, sigma2RM, vYwin, vYwinLag, vOls, dMi, init=<1;.5;0;0;2;30;1>, vPg=init[:2], vPst=init, dPg, dPst;

    if(cT<cWin) println("Catastrophic Error: rows(vY) < cWin!");
    s_obj=new Garch();
    for(i=0; i<=cT-cWin; ++i)
    {
        vYwin = vY[(i+1):(i+cWin-1)];
        vYwinLag = vY[i:(i+cWin-2)];
        olsc(vYwin, 1~vYwinLag, &vOls);
        dMi = (1~vYwin[cWin-2])*vOls;
        s_vE = vYwin-(1~vYwinLag)*vOls;
        if(MaxBFGS(LG, &vPg, &dPg, 0, 1))
        {
            vPg=init[:2];
            MaxBFGS(LG, &vPg, &dPg, 0, 1);
        }
        vPg[0]=max(0,vPg[0]);
        vPg[1]=max(0,vPg[1]);
        vPg[2]=min(max(0,vPg[2]),1);
        vPg[1]=min(1-vPg[2],vPg[1]);
        if(MaxBFGS(LST, &vPst, &dPst, 0, 1))
        {
            vPst=init;
            MaxBFGS(LST, &vPst, &dPst, 0, 1);
        }
        vPst[0]=max(0,vPst[0]);
        vPst[1]=max(0,vPst[1]);
        vPst[2]=min(max(-1,vPst[2]),1);
        vPst[3]=min(max(0,vPst[3]),1);
        vPst[1]=min(1-vPst[3],vPst[1]);
        vPst[4]=max(0,vPst[4]);
        vPst[5]=max(2,vPst[5]);
        vPst[6]=max(0,vPst[6]);
        sigma2RM = s_vE[0]^2;
        for(j=0; j<cWin-1; ++j) sigma2RM |= .06*s_vE[j]^2 + .94*sigma2RM[j];
        mVaR[i][0] = -(dMi + sqrt(sigma2RM[cWin-1])*quantnorm);
        mVaR[i][1] = -(dMi + sqrt(garch(vPg)[cWin-1])*quantnorm);
        mVaR[i][2] = -(dMi + sqrt(aparch(vPst[:4])[cWin-1])*s_obj.INVCDFTA(tau, log(vPst[6]), vPst[5]));
    }
    savemat("mVaR.mat", mVaR, 1);
}