% File mrsgarchestfor_con_all.m
% File to estimate Markov Regime-Switching GARCH
% Written by Juri Marcucci
% PhD student in Economics @
% Dept of Economics - UCSD
% Constrained optimization: fmincon
% The model is an M-state Markov Switching model with GARCH effects
% Main references:  Gray (1996, Journal of Financial Economics)
%                   Klaassen (2002, Empirical Economics)
% Author: Juri Marcucci
% Dept. of Economics - UCSD
%
% Revised version: January 2011
% This program has been tested on new versions of Matlab (2008a, 2009a, 2010a).

clc
clear all
format compact
warning off
seconds=clock;
%Loading data
[y]=textread('sugar.txt');

N = size(y,1);
T=2380; % In-Sample: 1/1/1988-9/28/2001
dataf = 100*(log(y(2:N))-log(y(1:N-1))); % data to be forecast
N = N - 1;
clear y;
%ERRORS
datafm = mean(dataf); %mean of all the sample for the forecast error
errors = dataf - datafm; % errors
data = dataf(1:T); % In-Sample data
% ERRORTYPE -> || 1 = NORMAL || 2 = t-DF || 4 = t 1df || 5 = GED 1df ||
% In the following loop, just put 'errortype = [j]', j=(1,2,4,5) to estimate only one model reducing
% the computational burden. You can leave it as it is, but it would take around 4/5 days to run everything.

for errortype = [1 2 4 5];
    data = dataf(1:T); % In-Sample data
    % Initialization of the parameters
    % Starting values
    % x0(1,1) = a01; x0(2,1) = a02; x0(3,1) = b01; x0(4,1) = b02; x0(5,1) = b11;
    % x0(6,1) = b12; x0(7,1) = b21; x0(8,1) = b22; x0(9,1) = p;  x0(10,1) = q;
    % x0(11,1) = \nu1; x0(12,1) = \nu2;
    x0=[1.5*datafm; %a01
            .5*datafm; %a02
            .0005; %b01
            .0002; %b02
            .15; %b11
            .25; %b12
            .6; %b21
            .7; %b22
            .9; %p
            .9  %q
            4.5; %\nu1
            4.5 %\nu2
            ];
    switch errortype
        case 1 %Normal
            x0=x0(1:10,1);
        case {2,3} % t & GED with 2DoF
            x0 = x0(1:12,1) ;
        case {4,5} % t and GED with 1DoF
            x0 = x0(1:11,1);
        otherwise
            error('Variable errortype must be integer in (1,5)')
    end
    original_param = zeros(size(x0,1),N-T+1);
    VF = [];
    Loglike = zeros(size(x0,1),N-T+1);
    Exitflag = zeros(size(x0,1),N-T+1);
    Parameters = zeros(size(x0,1),N-T+1,11);
    % Parameters for the conditional mean: a01, a02, i.e. x0(1 through 2)
    % Parameters for the conditional variance: b01, b02, b11, b12, b21, b22, i.e. x0(3 through 8)
    % Parameters for the P: p, q i.e. x0(9 through 10)
    % Parameters for the DoF: \nu_j, j=1,2, i.e. x0(11 through 12)
    % Print all results in a file
    modelname = ['MRSGARCH_'];
    switch errortype
        case 1
            distr_name = ['N'];
        case 2
            distr_name = ['t_2DoF'];
        case 4
            distr_name = ['t_1DoF'];
        case 5
            distr_name = ['GED_1DoF'];
        otherwise
    end
    % Mac/Linux
    filename=strcat('junk/',modelname,distr_name,'.txt')
    % Windows
    % filename=strcat('junk\',modelname,distr_name,'.txt')
    fid = fopen(filename,'w');
    tit1 = strcat('Model: \n',char(32),modelname,distr_name,' \n\n');
    fprintf(fid,'Results from optimization In Sample \n');
    fprintf(fid,tit1);

%%%%%%%%%%%%%%%% MAIN LOOP %%%%%%%%%%%%%%%%%%%%%%%%%%%
    index = 0; % time index for the rolling window
    constr = 0; % unconstrained paramters
    flag1=1 % 1 for Klaassen and 0 for Gray (default)
    max_forecast = 22; % max number of forecasts
    % Elapsed minutes from beginning
    minutes = etime(clock,seconds)/60;
    eval(['save junk/',strcat('res_',modelname,distr_name),' errortype constr']);

    for timeindex = T:N
        data = dataf(1+index:timeindex);
        sec=clock; % clock call for the main loop
        % Setting the lower bounds and upper bounds for optimization (once for all)
        % Using fmincon
        infinity = inf;
        hundred = 100;
        ten = 10;
        ptiny = 1e-6;
        tiny = 1e-8;
        ttiny = 1e-8;
        max_gamma = 171.6;
        ged_nu_tiny = 1/max_gamma;
        tdist_tiny = 1e-3; % New
   
        switch errortype
            case 1 % Normal case
                lowerBounds  =  [-hundred(ones(2,1)) ; tiny(ones(2,1)) ; ttiny(ones(2+2,1)); ptiny(ones(2,1))];
                upperBounds  =  [ hundred(ones(2,1)) ; ten(ones(2,1)) ;  ones(2+2,1) ; 1-ptiny(ones(2,1))];
            case {2} % t with 2 Dof
                lowerBounds  =  [-.5*hundred(ones(2,1)) ; tiny(ones(2,1)) ; tiny(ones(2+2,1)); ptiny(ones(2,1)) ; (2+tdist_tiny)+tiny(ones(2,1))]; % New
                upperBounds  =  [ .5*hundred(ones(2,1)) ; ten(ones(2,1)) ;  ones(2+2,1) ; 1-ptiny(ones(2,1)) ; 2*max_gamma-1+tiny(ones(2,1))];
            case 4 % t with 1 DoF
                lowerBounds  =  [-.5*hundred(ones(2,1)) ; tiny(ones(2,1)) ; zeros(2+2,1); ptiny(ones(2,1)) ; 2+tiny+tdist_tiny]; % New
                upperBounds  =  [ .5*hundred(ones(2,1)) ; ten(ones(2,1)) ;  ones(2+2,1) ; 1-ptiny(ones(2,1)) ; 2*max_gamma-1+tiny];
            case 5 % GED with 1 DoF
                lowerBounds  =  [-2*ten(ones(2,1)) ; tiny(ones(2,1)) ; tiny(ones(2+2,1)); ptiny(ones(2,1)) ; ged_nu_tiny*3+tiny];
                upperBounds  =  [ 2*ten(ones(2,1)) ; 2*(ones(2,1)) ; ones(2+2,1); 1-ptiny(ones(2,1)) ; max_gamma+tiny];
            otherwise
        end
        % Constraints for fmincon
        A = zeros(10,size(x0));
        A(1,5) = -1;
        A(2,6) = -1;
        A(3,7) = -1;
        A(4,8) = -1;
        A(5,9) = -1;
        A(6,9) = 1;
        A(7,10) = -1;
        A(8,10) = 1;
        A(9,[5 7]) = 1;
        A(10,[6 8]) = 1;
        b  =  [zeros(1,2) zeros(1,2) [0 1 0 1] ones(1,2)*[1  -  1e-6]];
        % For next call to the optimization procedure fmincon: Setting the options
        % Options for fmincon
        options  =  optimset('fmincon');
        if timeindex == T
            max_st_v = 99; % max # of starting values. Default=99
            switch errortype
                case 1 %N
                    values1 = [ x0 [abs(rand(size(x0,1)-2,500));.98*(ones(2,500))]];
                case 2 %t2DF
                    values1 = [ x0 [abs(rand(size(x0,1)-4,500));.98*(ones(2,500));2+6*abs(rand(2,500))]];
                case 4 %t1DF
                    values1 = [ x0 [abs(rand(size(x0,1)-3,500));.98*(ones(2,500));2+6*abs(rand(1,500))]];
                case 5 %G2DF
                    values1 = [ x0 [abs(rand(size(x0,1)-3,500));.98*(ones(2,500));abs(randn(1,500))]];
                otherwise
            end
            inndd1=(A(9:10,:)*values1>1);
            inndd2=sum(inndd1);
            inndd2(find(inndd2==2))=1;
            values1(:,(find(inndd2==1)))=[];
            starting_values = values1(:,1:max_st_v);
            clear inndd*;
            clear values1;
            clear values1;
            param0 = repmat(zeros(size(x0)),1,max_st_v);
            options  =  optimset(options , 'TolFun'      , 1e-7);
            options  =  optimset(options , 'TolX'        , 1e-7);
            options  =  optimset(options , 'TolCon'      , 1e-6);
            options  =  optimset(options , 'Display'     , 'iter');
            options  =  optimset(options , 'Diagnostics' , 'on');
            options  =  optimset(options , 'LargeScale'  , 'off');
            options  =  optimset(options , 'MaxIter'     , 2500);
            options  =  optimset(options , 'HessUpdate'     ,'bfgs');
            options  =  optimset(options , 'LevenbergMarquardt'     ,'on');
            options  =  optimset(options , 'LineSearchType'     ,'quadcubic');
            options  =  optimset(options , 'Jacobian'     ,'off');
            options  =  optimset(options , 'MeritFunction'     ,'singleobj');
            options  =  optimset(options , 'MaxFunEvals' , 2000);
            fprintf(fid,'Optimization for STARTING VALUES \n',index);
            LL=zeros(max_st_v,1);
            for j = 1:max_st_v
                fprintf(1,'Index is %-3.3f \n',j);
                fprintf(fid,'Index is %-3.3f \n',j);
                x0 = starting_values(:,j);
                fprintf(fid,'x0 \n');
                fprintf(fid,'%-3.5f \n',x0);
                LL(j,1)=swgarchlik(x0,data,errortype,flag1,constr);
                fprintf(fid,'LL is \n');
                fprintf(fid,'%-3.5f \n',LL(j,1));
                [parameters, LL(j,1)] =  fmincon('swgarchlik', x0 ,A, b, [],[],lowerBounds, upperBounds, [], options, data, errortype, flag1, constr);
                fprintf(fid,'LL is %-3.5f and param_x0 is \n',LL(j,1));
                fprintf(fid,'%-3.5f \n',parameters);
                param0(:,j) = parameters;
            end
            [minLL,minIndex] = min(LL,[],1)
            x0 = [];
            x0 = param0(:,minIndex)
            fprintf(fid,'Initial x0 \n');
            fprintf(fid,'%-3.5f \n',x0);
        end
        if timeindex == T
            options  =  optimset(options , 'TolFun'      , 1e-008);
            options  =  optimset(options , 'TolX'        , 1e-008);
        else            
            switch errortype
                case {1,5}
                    options  =  optimset(options , 'TolFun'      , 1e-008);
                    options  =  optimset(options , 'TolX'        , 1e-008);
                otherwise
                    options  =  optimset(options , 'TolFun'      , 1e-007);
                    options  =  optimset(options , 'TolX'        , 1e-007);
            end
        end
        options  =  optimset(options , 'TolCon'        , 1e-008);
        options  =  optimset(options , 'TolPCG'        , 1e-2);
        options  =  optimset(options , 'Display'     , 'iter');
        options  =  optimset(options , 'Diagnostics' , 'on');
        options  =  optimset(options , 'LargeScale'  , 'off');
        options  =  optimset(options , 'MaxIter'     , 2500);
        options  =  optimset(options , 'HessUpdate'     ,'bfgs');
        options  =  optimset(options , 'LevenbergMarquardt'     ,'off');
        options  =  optimset(options , 'LineSearchType'     ,'cubicpoly');
        options  =  optimset(options , 'Jacobian'     ,'off');
        options  =  optimset(options , 'MeritFunction'     ,'multiobj');
        options  =  optimset(options , 'MaxFunEvals' , 5000);
        
        fprintf(1,'\n\n Index is %-3.0f \n',index);
        fprintf(1,'\n ERRORTYPE is %-3.0f \n',errortype);
        fprintf(1,'\n Minutes from latest optimization are %-3.2f \n',minutes);
        fprintf(1,'Time from latest optimization is %-15s \n',datestr(clock,0));
        
        [parameters, LLF, EXITFLAG, OUTPUT, Lambda, GRAD, HESSIAN] =  fmincon('swgarchlik', x0 ,A, b, [],[],lowerBounds, upperBounds, [], options, data, errortype, flag1, constr);
   
        if EXITFLAG<=0
            EXITFLAG
            fprintf(1,'Convergence has been not successful! \n')
        end
        % Add check for parameter constraints if fminunc is used
        % Printing
        fprintf(fid,'timeindex \t\t LogLike \t\t Exitflag(<=0_if_unsuccessful) \n');
        fprintf(fid,'%-3.1f \t\t %-5.5f \t\t %-3.3f \n',[timeindex,LLF,EXITFLAG]);
        [LLF,likelihoods,e,p1,p1t,p1sm,h1,h2] = swgarchlik(parameters,data,errortype,flag1,constr);
        LLF=-LLF;
        grad = gradt('swgarchlik',parameters,data,errortype,flag1,constr);
        hess = hessian('swgarchlik',parameters,data,errortype,flag1,constr);
        % asymptotic standard errors
        stder_hess = sqrt(diag((hess)^(-1))); % from calculated hessian
        tstat_hess = parameters./stder_hess;

%%%%%%%%%%%%%%%%%%%%
        % ADDED THE FOLLOWING 2 LINES FOR STD ERR WITH HESSIAN FROM
        % OPTIMIZATION. they were disappeared and I don't know why
        stder_HESS = sqrt(diag((HESSIAN)^(-1)));
        tstat_HESS = parameters./stder_HESS;
        % robust standard errors
        h = min(abs(parameters/2),max(parameters,1e-2))*eps^(1/3);
        hplus=parameters+h;
        hminus=parameters-h;
        likelihoodsplus=zeros(T,length(parameters));
        likelihoodsminus=zeros(T,length(parameters));
        for i=1:length(parameters)
            hparameters=parameters;
            hparameters(i)=hplus(i);
            [holder1,indivlike] = swgarchlik(hparameters,data,errortype,flag1,constr);
            likelihoodsplus(:,i)=indivlike;
        end
        for i=1:length(parameters)
            hparameters=parameters;
            hparameters(i)=hminus(i);
            [holder1,indivlike] = swgarchlik(hparameters,data,errortype,flag1,constr);
            likelihoodsminus(:,i)=indivlike;
        end
        clear holder*
        scores=(likelihoodsplus-likelihoodsminus)./(2*repmat(h',T,1));
        scores=scores-repmat(mean(scores),T,1);
        B=(1/T)*scores'*scores;
        OPG=scores'*scores;
        A=(1/T)*hess;
        SE_rob_hess=sqrt((1/T)*diag((A^(-1))*B*(A^(-1)))); % hessian from calculation
        %%%%%%%%%%%%%%%%%%%%%%%%%
        % ADDED THE FOLLOWING 3 LINES FOR SE_rob_HESS
        A2=(1/T)*HESSIAN;
        SE_rob_HESS=sqrt((1/T)*diag((A2^(-1))*B*(A2^(-1)))); % hessian from OPTIMIZATION
        tstat_HESS_rob = parameters./SE_rob_HESS; % hessian from calculation
        clear likelihoods*;
        clear indivlike*;
        clear scores;
        tstat_hess_rob = parameters./SE_rob_hess; % hessian from calculation
        stderrors = sqrt(diag((OPG^(-1)))); %Inverse of OPG (outer product of gradients)
        tstats = parameters./stderrors;
        % Printing everything
        % COMMENTED THE NEXT 5 LINES. THEY WOULD WRONGLY SET THE GRAD, STD
        % ERR AND TSTAT TO [] BEFORE STORING THEM!
        %GRAD=[];
        %stder_HESS=[];
        %SE_rob_HESS=[];
        %tstat_HESS=[];
        %tstat_HESS_rob=[];
        fprintf(fid,'Param \t\t GRAD1(optim) \t\t Grad2(gradt) \t\t stdH \t\t stdh \t\t stdRH \t\t stdRh \t\t stderr \t\t tH \t\t th  \t\t tHR \t\t thR \t\t tstats \n');
        fprintf(fid,'%-3.6f \t\t %-3.6f \t\t %-3.3f \t\t %-3.3f \t\t %-3.3f \t\t %-3.3f \t\t %-3.3f \t\t %-3.3f \t\t %-3.3f \t\t %-3.3f \t\t %-3.3f \t\t %-3.3f \t\t %-3.3f \n',[parameters(:),GRAD(:),grad(:),stder_HESS(:),stder_hess(:),SE_rob_HESS(:),SE_rob_hess(:),stderrors(:),tstat_HESS(:),tstat_hess(:),tstat_HESS_rob(:),tstat_hess_rob(:),tstats(:)]');

        % Cheching if parameters are on the boundaries
        original_param(:,index+1) = parameters(:);
        index_lb = find(parameters == lowerBounds);
        index_ub = find(parameters == upperBounds);
        parameters(index_lb) = parameters(index_lb).*rand(size(index_lb));
        parameters(index_ub) = parameters(index_ub).*rand(size(index_ub));

        %%%%%%%%%%%%%%%% STORING RESULTS %%%%%%%%%%%%%%%%%%%%%%%%%%
      
        Loglike(:,index+1) = LLF;
        Exitflag(:,index+1) = EXITFLAG;
   
        Parameters(:,index+1,1) = parameters(:);
        Parameters(:,index+1,2) = GRAD(:);
        Parameters(:,index+1,3) = grad(:);
        Parameters(:,index+1,4) = stder_HESS(:);
        Parameters(:,index+1,5) = stder_hess(:);
        Parameters(:,index+1,6) = SE_rob_HESS(:);
        Parameters(:,index+1,7) = SE_rob_hess(:);
        Parameters(:,index+1,8) = tstat_HESS(:);
        Parameters(:,index+1,9) = tstat_hess(:);
        Parameters(:,index+1,10) = tstat_HESS_rob(:);
        Parameters(:,index+1,11) = tstat_hess_rob(:);
   
        minutes=etime(clock,sec)/60;
        % Saving
        eval(['save junk/',strcat('res_',modelname,distr_name),' VF Loglike Exitflag Parameters minutes index -append']);

        %%%%%%%%%%% SETTING THE NEW STARTING POINT %%%%%%%%%%%%%%%%%
        x0 = parameters; %starting value = previous optimal value. Otherwise
        %comment and you will start from same starting value every time.
        clear data;
        index = index + 1 % index for the rolling window
        % Elapsed minutes since last optimization
        minutes = etime(clock,sec)/60;
    end
    % Total time
    minutes = etime(clock,seconds)/60;
    hours = etime(clock,seconds)/3600;
    fprintf(1,'Minutes \t\t Hours \n')
    fprintf(fid,'Minutes \t\t Hours \n')
    fprintf(1,'%-3.3f \t\t %-3.3f \n',[minutes hours]')
    fprintf(fid,'%-3.3f \t\t %-3.3f \n',[minutes hours]')

    status = fclose(fid);

    modelsavename = ['mrsgarch_'];
    switch errortype
        case 1
            distr_name = ['n'];
        case 2
            distr_name = ['t2'];
        case 4
            distr_name = ['t'];
        case 5
            distr_name = ['g'];
        otherwise
    end
    eval(['save ',strcat(modelsavename,distr_name)]);
end

上面是Juri Marcucci的MRS-GARCH的程序

好像发现了一点头绪，这个var_forec_mrsg_var.m里的T值也要修改