
function [parameters, LLF, log_likelihood, ht, smooth_t, filter_t, stderrors]=MLE_RS_Uni_GARCH_fit(data)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Estimate the Univariate GARCH (1,1) model under regime switching framework, follows Gray (1996)
%% since under regime switching, GARCH(1,1) is enough in empirical analysis
% Written By Weijun Xu, 
% 01/08/2008
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Let's first set some initial values, alpha(st)+beta(st)<1 and alpha(st),
%% beta(st)>0, under each regime. 0<tp11,tp00<1
% we add a little prior informaiton on the constant term but not on the alpha,
% and beta, since it is hard to say which regime has higher or lower alpha and beta
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T=length(data);

alpha=0.05;
beta=0.92;
omeg1=(1-(alpha+beta))*1.5*var(data);    %Use variance targeting form to its initial value at regime 1 (with higher variance)
omeg0=(1-(alpha+beta))*0.5*var(data);    %Use variance targeting form to its initial value at regime 0 (with lower variance)
tp11=0.9;
tp00=0.9;

para_start=[omeg1,alpha,beta, omeg0,alpha,beta, tp11, tp00]'; 

num_para=length(para_start);
Aineq=[0,1,1,zeros(1,5); zeros(1,3),0,1,1,zeros(1,2)];
Cineq=[(1-1e-3)*ones(2,1)];
LB=[(1e-3)*ones(num_para,1)];
UB=[inf*ones(num_para-2,1); (1-1e-6)*ones(2,1)];

options=optimset('fmincon');
options=optimset(options , 'Display','iter','Diagnostics','off','LargeScale','off','MaxIter',2000, 'MaxFunEvals', 1e+6);

[parameters]=...
   fmincon(@(parameters)MLE_RS_Uni_GARCH_likelihood(parameters,data),para_start,Aineq,Cineq,[],[],LB,UB,[],options);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Calculate Kim's Smooth probability: Smooth_t(t)=filter_t(t).*{[Smooth_t(t+1)./Conditional_t(t+1)]*P'};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[LLF,log_likelihood, ht, filter_t,conditional_t]=MLE_RS_Uni_GARCH_likelihood(parameters, data);
LLF=-LLF;

P=[parameters(end-1), 1-parameters(end-1); 1-parameters(end), parameters(end)];
smooth_t=zeros(T,2);
smooth_t(end,:)=filter_t(end,:);
for i=T-1:-1:1
    smooth_t(i,:)=filter_t(i,:).*((smooth_t(i+1,:)./conditional_t(i+1,:))*P');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Calculate the standard errors, and return conditional varianec {ht}, high volatility variance {ht_1}, low volatility variance {ht_0} 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the Hessian

if nargout>5
    fprintf('\n Calculating the standard errors of MRS-Uni-GARCH model.... \n')
    hessian=hessian_2sided('MLE_RS_Uni_GARCH_likelihood',parameters, data);
    h=max(abs(parameters/2),1e-2)*eps^(1/3);
    hplus=parameters+h;
    hminus=parameters-h;
    likelihoodsplus=zeros(T,length(parameters));
    likelihoodsminus=zeros(T,length(parameters));
    % Calculate the Gradients
    for i=1:length(parameters)
        hparameters=parameters;
        hparameters(i)=hplus(i);
        [HOLDER, indivlike] =MLE_RS_Uni_GARCH_likelihood(hparameters, data);
        likelihoodsplus(:,i)=indivlike;
    end
    for i=1:length(parameters)
        hparameters=parameters;
        hparameters(i)=hminus(i);
        [HOLDER, indivlike] =MLE_RS_Uni_GARCH_likelihood(hparameters, data);
        likelihoodsminus(:,i)=indivlike;
    end
    scores=(likelihoodsplus-likelihoodsminus)./(2*repmat(h',T,1));
    B=scores'*scores;
    stderrors=inv(hessian)*B*inv(hessian);
    stderrors=sqrt(diag(stderrors));
end

% Function End