% DO_OUT_SAMPLE.M       Performs the out-of-sample comparison of the 
%   one-step-ahead Expected Shrotfall (ES) forecasts (1%, 5% and 10% ES) 
%   obtained from three competing models: (1) APD-GARCH(1,1) model from 
%   Equations (10)-(11) in Komunjer (2006), (2) asymmetric Student-t-GARCH(1,1) 
%   model from Equation (12) in Komunjer (2006), and (3) GPD-EGARCH(1,1) model
%   from Equations (13)-(14) in Komunjer (2006).
%
% References:
%   [1] Komunjer, I. (2006): "Asymmetric Power Distribution: Theory and 
%   Applications to Risk Measurement", Journal of Applied Econometrics (forthcoming)
%
% Author: Ivana Komunjer - UCSD (komunjer@ucsd.edu)
% Version: 2.1        Date: 09/11/2006

global GARCH_SERIES RES_SERIES;
datatype = ['S&P500_1990_2002.txt'; 'NASDAQ_1990_2002.txt'; ...
        'MSFT_1990_2002.txt  ';'EXUKUS_1990_2002.txt'];
datatypeout = ['S&P500_2002_2006.txt'; 'NASDAQ_2002_2006.txt'; ...
        'MSFT_2002_2006.txt  ';'EXUKUS_2002_2006.txt'];
type = cellstr(datatype); typeout = cellstr(datatypeout);
Ntype = 1;
while Ntype <= size(type,1)
    % IN-SAMPLE analysis 
    timeseries = char(type(Ntype))
    [DATE, PRICES] = textread(timeseries,'%u%n', 'headerlines', 4);
    raw_xlsdate = DATE(2:size(DATE,1));            % date in xls format - convention 1900
    raw_series = 100*(log(PRICES(2:size(PRICES,1)))-log(PRICES(1:size(PRICES,1)-1)));
    indSplit = logical(ones(size(raw_xlsdate,1),1));
    if Ntype == 3          % removes split dates for MSFT
        indSplit = (raw_xlsdate ~= 32975)&(raw_xlsdate ~= 33416)&(raw_xlsdate ~= 33770)&(raw_xlsdate ~= 34477)...
            &(raw_xlsdate ~= 35408)&(raw_xlsdate ~= 35849)&(raw_xlsdate ~= 36248);
    end
    raw_series = raw_series(indSplit);
    raw_xlsdate = raw_xlsdate(indSplit);
    isnozero = (isfinite(raw_series));        
    xlsdate = raw_xlsdate(isnozero);
    GARCH_SERIES = raw_series(isnozero);
    matdate = x2mdate(xlsdate, 0);      % converts xls format into Matlab format date
    Yt = GARCH_SERIES;
    % OUT-OF-SAMPLE analysis
    timeseriesout = char(typeout(Ntype));
    [DATEout, PRICESout] = textread(timeseriesout,'%u%n', 'headerlines', 2);
    raw_xlsdateout = DATEout(2:size(DATEout,1));            % date in xls format - convention 1900
    raw_seriesout = 100*(log(PRICESout(2:size(PRICESout,1)))-log(PRICESout(1:size(PRICESout,1)-1)));
    indSplitout = logical(ones(size(raw_xlsdateout,1),1));
    raw_seriesout = raw_seriesout(indSplitout);
    raw_xlsdateout = raw_xlsdateout(indSplitout);
    isnozeroout = (isfinite(raw_seriesout));         
    xlsdateout = raw_xlsdateout(isnozeroout);
    GARCHout = raw_seriesout(isnozeroout);
    matdateout = x2mdate(xlsdateout, 0);      % converts xls format into Matlab format date
    Ytout = GARCHout;
    plot(matdateout, Ytout);
    axis([matdateout(1) matdateout(size(matdateout,1)) -6 6]);
    dateaxis('x',12);

    % APD MAX-LIK
    alpha10 = 0.8; if Ntype == 4 alpha10 = 0.7; end % for the BP/USD rate series!!!
    alpha20 = 0.05; alpha00 = (1-alpha10-alpha20)*(std(Yt)^2); 
    mu0 = mean(Yt);
    alpha0 = 0.5; lda0 = 2;
    par0 = [alpha00; alpha10; alpha20; mu0; alpha0; lda0];
    [parML,loglik,exitflag,output,lambda,grad,hessian] = MLAPDGARCH(par0);
    % computing the alphabar-expected shortfall and alpha-bar quantile
    ESepsall = []; quantepsall = [];
    alphabar = [0.01; 0.05; 0.10];
    for i=1:1:size(alphabar,1)
        aAQ = parML(5); 
        ldaAQ = parML(6);
        deltaAQ = [2*(aAQ^ldaAQ)*((1-aAQ)^ldaAQ)]/[(aAQ^ldaAQ)+((1-aAQ)^ldaAQ)];
        muAQ = (gamma(2/ldaAQ))/[(gamma(3/ldaAQ))*(gamma(1/ldaAQ))*(1-3*aAQ+3*aAQ^2)-(gamma(2/ldaAQ))^2*(1-2*aAQ)^2]^0.5;
        thetaAQ = -1*(1-2*aAQ)*muAQ;
        phiAQ = gamma(1/ldaAQ)*(deltaAQ^(1/ldaAQ))*muAQ/gamma(2/ldaAQ);
        quanteps = thetaAQ + phiAQ*AlphaQuantInv(aAQ,ldaAQ,alphabar(i));
        ESeps = phiAQ*APDES(alphabar(i),aAQ,ldaAQ);
        ESepsall = [ESepsall; ESeps]; quantepsall = [quantepsall; quanteps];
    end
    % forecasting APD volatility
    T = size(Yt,1);
    % reconstructing the volatility at time T
    sig2tML = parML(1)/(1-parML(2)-parML(3))*(1:1:T)';
    for i = 2:1:T
        sig2tML(i) = parML(1) + parML(3)*sig2tML(i-1) + parML(2)*((Yt(i-1)-parML(4))^2);
    end
    R = size(Ytout,1);
    sig2tout = (parML(1) + parML(3)*sig2tML(T) + parML(2)*(Yt(T)-parML(4))^2)*(1:1:R)';
    for j = 2:1:R
        sig2tout(j) = parML(1) + parML(3)*sig2tout(j-1) + parML(2)*((Ytout(j-1)-parML(4))^2);
    end
    % computing the residuals
    epstout = (Ytout - parML(4))./(sig2tout.^0.5);
    % computing the ES residual statistic
    Reps = ones(R,1)*quantepsall' - epstout*ones(1,3) - ones(R,1)*ESepsall';
    IndRviol = (epstout*ones(1,3) <= ones(R,1)*quantepsall');
    % plots residuals R and violations for 5%
    subplot(3,1,1);
    plot(matdateout, Reps(:,2), 'b', matdateout(IndRviol(:,2)),Reps(IndRviol(:,2),2),'k.');
    axis([matdateout(1) matdateout(size(matdateout,1)) -7 2]);
    dateaxis('x',12);
    meanReps = []; stdReps = []; meanquanteps = [];
    for ii = 1:1:3
        Rfinal = Reps(IndRviol(:,ii),ii); 
        meanReps = [meanReps; mean(Rfinal)]; meanquanteps = [meanquanteps; mean(IndRviol(:,ii))];
        stdReps = [stdReps; std(Rfinal)];
    end
    % saves the APD results
    APDstat = [meanReps stdReps meanReps./stdReps ESepsall meanquanteps quantepsall];
    
    % Skew-t MAX-LIK
    alpha10 = 0.8; if Ntype == 4 alpha10 = 0.7; end % for the BP/USD rate series!!!
    alpha20 = 0.05; alpha00 = (1-alpha10-alpha20)*(std(Yt)^2); 
    mu0 = mean(Yt);
    xi0 = 0; nu0 = 5;
    par0 = [alpha00; alpha10; alpha20; mu0; xi0; nu0];
    [parML,loglik,exitflag,output,lambda,grad,hessian] = MLSkewTGARCH(par0);
    % computing the alphabar-expected shortfall and alpha-bar quantile
    ESepsall = []; quantepsall = [];
    alphabar = [0.01; 0.05; 0.10];
    for i=1:1:size(alphabar,1)
        nuT = parML(6); 
        ldaT = ((exp(parML(5)))^2-1)/((exp(parML(5)))^2+1);
        quanteps = skewtdis_inv(alphabar(i), nuT, ldaT);
        % computing ES (my computations based on Giot and Laurent 2004)
        ESeps = SkewTES(alphabar(i),exp(parML(5)),nuT);
        ESepsall = [ESepsall; ESeps]; quantepsall = [quantepsall; quanteps];
    end
    % forecasting skewT volatility
    T = size(Yt,1);
    % reconstructing the volatility at time T
    sig2tML = parML(1)/(1-parML(2)-parML(3))*(1:1:T)';
    for i = 2:1:T
        sig2tML(i) = parML(1) + parML(3)*sig2tML(i-1) + parML(2)*((Yt(i-1)-parML(4))^2);
    end
    R = size(Ytout,1);
    sig2tout = (parML(1) + parML(3)*sig2tML(T) + parML(2)*(Yt(T)-parML(4))^2)*(1:1:R)';
    for j = 2:1:R
        sig2tout(j) = parML(1) + parML(3)*sig2tout(j-1) + parML(2)*((Ytout(j-1)-parML(4))^2);
    end
    % computing the residuals
    epstout = (Ytout - parML(4))./(sig2tout.^0.5);
    % computing the ES residual statistic
    Reps = ones(R,1)*quantepsall' - epstout*ones(1,3) - ones(R,1)*ESepsall';
    IndRviol = (epstout*ones(1,3) <= ones(R,1)*quantepsall');
    % plots residuals R and violations for 5%
    subplot(3,1,2);
    plot(matdateout, Reps(:,2), 'b', matdateout(IndRviol(:,2)),Reps(IndRviol(:,2),2),'k.');
    axis([matdateout(1) matdateout(size(matdateout,1)) -7 2]);
    dateaxis('x',12);
    meanReps = []; stdReps = []; meanquanteps = [];
    for ii = 1:1:3
        Rfinal = Reps(IndRviol(:,ii),ii); 
        meanReps = [meanReps; mean(Rfinal)]; meanquanteps = [meanquanteps; mean(IndRviol(:,ii))];
        stdReps = [stdReps; std(Rfinal)];
    end
    % saves the SkewT results
    skewTstat = [meanReps stdReps meanReps./stdReps ESepsall meanquanteps quantepsall];

    % asymmetric GPD MAX-LIK
    alpha10 = 0.8; % alpha10 = 0.7; % for the BP/USD rate series!!!
    alpha20 = 0.05; 
    alpha30 = 0.05;
    alpha00 = (1-alpha10)*2*log(std(Yt)); 
    mu0 = mean(Yt);
    alpha0 = 0.5; lda0 = 2;
    par0 = [alpha00; alpha10; alpha20; alpha30; mu0; alpha0; lda0];
    [parML,loglik,exitflag,output,lambda,grad,hessian] = MLGPDEGARCH(par0);
    % computing the alphabar-expected shortfall and alpha-bar quantile
    ESepsall = []; quantepsall = [];
    alphabar = [0.01; 0.05; 0.10];
    for i=1:1:size(alphabar,1)
        aAQ = 0.5; 
        ldaAQ = parML(7);
        deltaAQ = [2*(aAQ^ldaAQ)*((1-aAQ)^ldaAQ)]/[(aAQ^ldaAQ)+((1-aAQ)^ldaAQ)];
        muAQ = (gamma(2/ldaAQ))/[(gamma(3/ldaAQ))*(gamma(1/ldaAQ))*(1-3*aAQ+3*aAQ^2)-(gamma(2/ldaAQ))^2*(1-2*aAQ)^2]^0.5;
        thetaAQ = -1*(1-2*aAQ)*muAQ;
        phiAQ = gamma(1/ldaAQ)*(deltaAQ^(1/ldaAQ))*muAQ/gamma(2/ldaAQ);
        quanteps = thetaAQ + phiAQ*AlphaQuantInv(aAQ,ldaAQ,alphabar(i));
        ESeps = phiAQ*APDES(alphabar(i),aAQ,ldaAQ);
        ESepsall = [ESepsall; ESeps]; quantepsall = [quantepsall; quanteps];
    end
    % forecasting GPD volatility
    T = size(Yt,1);
    % reconstructing the volatility at time T
    epsMLt = GARCH_SERIES - parML(4);
    sig2tML = (std(Yt)^2)*(1:1:T)';
    for i = 2:1:T
        sig2tML(i) = exp(parML(1) + parML(2)*log(sig2tML(i-1)) + parML(3)*(epsMLt(i-1)) + parML(4)*abs(epsMLt(i-1)));
    end   
    R = size(Ytout,1);
    sig2tout = exp(parML(1) + parML(2)*log(sig2tML(T)) + parML(3)*(epsMLt(T)) + parML(4)*abs(epsMLt(T)))*(1:1:R)';
    epsMLtout = Ytout - parML(4);
    for j = 2:1:R
        sig2tout(j) = exp(parML(1) + parML(2)*log(sig2tout(j-1)) + parML(3)*(epsMLtout(j-1)) + parML(4)*abs(epsMLtout(j-1)));
    end
    % computing the residuals
    epstout = (Ytout - parML(4))./(sig2tout.^0.5);
    % computing the ES residual statistic
    Reps = ones(R,1)*quantepsall' - epstout*ones(1,3) - ones(R,1)*ESepsall';
    IndRviol = (epstout*ones(1,3) <= ones(R,1)*quantepsall');
    % plots residuals R and violations for 5%
    subplot(3,1,3);
    plot(matdateout, Reps(:,2), 'b', matdateout(IndRviol(:,2)),Reps(IndRviol(:,2),2),'k.');
    axis([matdateout(1) matdateout(size(matdateout,1)) -7 2]);
    dateaxis('x',12);
    meanReps = []; stdReps = []; meanquanteps = [];
    for ii = 1:1:3
        Rfinal = Reps(IndRviol(:,ii),ii); 
        meanReps = [meanReps; mean(Rfinal)]; meanquanteps = [meanquanteps; mean(IndRviol(:,ii))];
        stdReps = [stdReps; std(Rfinal)];
    end
    GPDstat = [meanReps stdReps meanReps./stdReps ESepsall meanquanteps quantepsall];
    % prints out the results
    	filename = ['OutSampleESforecasts' int2str(Ntype) '.txt'];
	fid = fopen(filename,'a');
	fprintf(fid,'%s  %s  %s\n','================',timeseries,'================');
	fprintf(fid,'\n');
    fprintf(fid,'\n');
	fprintf(fid,'%s\n','APD forecasts');
	fprintf(fid,'\n');
	fprintf(fid,'%s\n','alpha_bar  mean(R) std(R) t-stat(R) ES(eps) mean(IndVaR) quant(eps):');
	fprintf(fid,'%6.3f  %6.3f  %6.3f  %6.3f  %6.3f  %6.3f  %6.3f\n',alphabar(1), APDstat(1,1), APDstat(1,2), APDstat(1,3), APDstat(1,4), APDstat(1,5), APDstat(1,6));
	fprintf(fid,'%6.3f  %6.3f  %6.3f  %6.3f  %6.3f  %6.3f  %6.3f\n',alphabar(2), APDstat(2,1), APDstat(2,2), APDstat(2,3), APDstat(2,4), APDstat(2,5), APDstat(2,6));
	fprintf(fid,'%6.3f  %6.3f  %6.3f  %6.3f  %6.3f  %6.3f  %6.3f\n',alphabar(3), APDstat(3,1), APDstat(3,2), APDstat(3,3), APDstat(3,4), APDstat(3,5), APDstat(3,6));
    fprintf(fid,'\n');
    fprintf(fid,'%s\n','Skew-T forecasts');
	fprintf(fid,'\n');
	fprintf(fid,'%s\n','alpha_bar  mean(R) std(R) t-stat(R) ES(eps) mean(IndVaR) quant(eps):');
	fprintf(fid,'%6.3f  %6.3f  %6.3f  %6.3f  %6.3f  %6.3f  %6.3f\n',alphabar(1), skewTstat(1,1), skewTstat(1,2), skewTstat(1,3), skewTstat(1,4), skewTstat(1,5), skewTstat(1,6));
	fprintf(fid,'%6.3f  %6.3f  %6.3f  %6.3f  %6.3f  %6.3f  %6.3f\n',alphabar(2), skewTstat(2,1), skewTstat(2,2), skewTstat(2,3), skewTstat(2,4), skewTstat(2,5), skewTstat(2,6));
	fprintf(fid,'%6.3f  %6.3f  %6.3f  %6.3f  %6.3f  %6.3f  %6.3f\n',alphabar(3), skewTstat(3,1), skewTstat(3,2), skewTstat(3,3), skewTstat(3,4), skewTstat(3,5), skewTstat(3,6));
	fprintf(fid,'\n');
    fprintf(fid,'\n');
	fprintf(fid,'%s\n','GPD forecasts');
	fprintf(fid,'\n');
	fprintf(fid,'%s\n','alpha_bar  mean(R) std(R) t-stat(R) ES(eps) mean(IndVaR) quant(eps):');
	fprintf(fid,'%6.3f  %6.3f  %6.3f  %6.3f  %6.3f  %6.3f  %6.3f\n',alphabar(1), GPDstat(1,1), GPDstat(1,2), GPDstat(1,3), GPDstat(1,4), GPDstat(1,5), GPDstat(1,6));
	fprintf(fid,'%6.3f  %6.3f  %6.3f  %6.3f  %6.3f  %6.3f  %6.3f\n',alphabar(2), GPDstat(2,1), GPDstat(2,2), GPDstat(2,3), GPDstat(2,4), GPDstat(2,5), GPDstat(2,6));
	fprintf(fid,'%6.3f  %6.3f  %6.3f  %6.3f  %6.3f  %6.3f  %6.3f\n',alphabar(3), GPDstat(3,1), GPDstat(3,2), GPDstat(3,3), GPDstat(3,4), GPDstat(3,5), GPDstat(3,6));
	fprintf(fid,'\n');
	fprintf(fid,'\n');
	fclose(fid);
    % next series
    Ntype = Ntype + 1;
end
