% DO_APD_MLE.M        Maximum Likelihood Estimation of the parameters of 
%       a GARCH(1,1) process with APD innovations:
%       Y[t] = mu + sig[t]*eta[t] where
%       sig[t]^2 = alpha0 + alpha1*sig[t-1]^2 + alpha2*(Y[t-1]-mu)^2 and
%       eta[t] is IID APD distributed with mean zero and unit variance (the pdf
%       of eta[t] is given by Equation (3) 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;
% Initialisation
datatype = ['S&P500_1990_2002.txt'; 'NASDAQ_1990_2002.txt'; ...
         'MSFT_1990_2002.txt  ';  'EXUKUS_1990_2002.txt'];
type = cellstr(datatype);
Ntype = 1;
while Ntype <= size(type,1)
    resmG = [];
    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;
    
    % Maximum Likelihood Estimation
    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);
    T = size(Yt,1);
    covMatrixML = inv(hessian)/T;
    stdML = diag(covMatrixML,0).^0.5;
    % computing the alphabar-expected shortfall and its standard deviation
    % by DELTA-METHOD
    alphabar = (0.005:0.005:0.995)';
    meanESYtall = ESr(Ntype,alphabar,parML);
    Nsim = 10;
    simpar = parML*ones(1,Nsim) + stdML*(2*(1:1:Nsim)/Nsim-1);
    indOK = ((simpar(2,:) + simpar(3,:))<1)&(simpar(1,:)> 0)&(simpar(2,:)> 0)&(simpar(6,:)> 0)&(simpar(5,:)> 0)&(simpar(5,:)< 1);
    ESYtgrad_all = [];
    for indpar = 1:1:6
        xpar = []; ESYtall = [];
        for iisim = 1:1:Nsim
            if indOK(iisim) == 1
                par = parML; par(indpar) = simpar(indpar,iisim); 
                ESYtpar = ESr(Ntype,alphabar,par);
                ESYtall = [ESYtpar, ESYtall];
                xpar = [xpar; simpar(indpar,iisim)];
            end
        end
        ESYtgrad_par_all = [];
        for indalpha = 1:1:size(alphabar,1)
            ESYtcol = ESYtall(indalpha,:); ESYtcol = ESYtcol';
            [xparsort, indsort] = sort(xpar);
            ESYtcolsort = ESYtcol(indsort);
            ESYtcol_grad = central_diff(ESYtcolsort, xparsort );
            indMLpar = (xparsort == parML(indpar));
            ESYtgrad_par = ESYtcol_grad(indMLpar);
            ESYtgrad_par_all = [ESYtgrad_par_all; ESYtgrad_par];
        end
        ESYtgrad_all = [ESYtgrad_all, ESYtgrad_par_all];
    end
    ESYtcovall = diag(ESYtgrad_all*covMatrixML*(ESYtgrad_all'));
    plot(alphabar, meanESYtall,'r',alphabar, meanESYtall-1.96*ESYtcovall.^0.5,'b--',alphabar, meanESYtall+1.96*ESYtcovall.^0.5,'b--');
    if Ntype == 1 SPESYtall = [meanESYtall ESYtcovall.^0.5]; elseif Ntype == 2 NASDESYtall = [meanESYtall ESYtcovall.^0.5];
    elseif Ntype == 3 MSFTESYtall = [meanESYtall ESYtcovall.^0.5]; elseif Ntype == 4 BPUSDESYtall = [meanESYtall ESYtcovall.^0.5]; end
    Ntype = Ntype + 1;
end
% saves the output
save ESstdDELTA SPESYtall NASDESYtall MSFTESYtall BPUSDESYtall;
% produces Figure 8 in the paper
plot(alphabar,SPESYtall(:,1)-1.96*SPESYtall(:,2),'k-',...
    alphabar,NASDESYtall(:,1)-1.96*NASDESYtall(:,2),'k--',...
    alphabar,MSFTESYtall(:,1)-1.96*MSFTESYtall(:,2),'k:',...
    alphabar,BPUSDESYtall(:,1)-1.96*BPUSDESYtall(:,2),'k.-',...
    alphabar,SPESYtall(:,1)+1.96*SPESYtall(:,2),'k-',...
    alphabar,NASDESYtall(:,1)+1.96*NASDESYtall(:,2),'k--',...
    alphabar,MSFTESYtall(:,1)+1.96*MSFTESYtall(:,2),'k:',...
    alphabar,BPUSDESYtall(:,1)+1.96*BPUSDESYtall(:,2),'k.-');
legend('S&P500', 'NASDAQ','Microsoft','BP/USD')
axis([0 1 0 3])
title('Delta Method 95%-Confidence Intervals for Expected Shortfalls on 3/16/1993')

% prints stuff on the screen
indout = [2;10;20;40;100];
[alphabar(indout) SPESYtall(indout,:) NASDESYtall(indout,:) MSFTESYtall(indout,:) BPUSDESYtall(indout,:)]
