ms_ar_fit.m

Autoregressive Markov Switching Model函数用于评估、仿真及预测自回归的马尔可夫转换模型。可以选择用于模型估计的分布函数。用于研究时间序列结构性变化

% Function for estimation of a Autoregressive Markov Switching model with k
% states and p lags (MS(k)-AR(p))
%
%   The models swicthes in the whole mean equation including AR coefficients and constant
%
%
%   Input:  x - Time series to be modelled
%           ar - Order of autoregressive component
%           k - Number of states
%           advOpt  - A structure with advanced options. The field advOpt.distrib is the distribution assumption ('Normal' or 't',
%                     default='Normal'). the field advOpt.std_method is the
%                     method for calculating the standard error of the
%                     fitted coefficients (1 for using second partial derivatives, 2 for using
%                     outer product matrix, 3 for using White's covariance
%                     matrix and 4 for using newey and West Matrix). More details about each can be found in the pdf file
%                     The default value is 1.
%                     The field useMex is for wheter to use mex version of filter in the likelihood calculation or
%                     not (default is no use of mex function, that is
%                     useMex=0).
%
%   Output: Spec_Output - A structure with followning fields:
%
%               LL - Log likelihood of fitted model
%
%               Probs - States probabilities over time (each collum represents
%               each state, ascending order).
%
%               Coeff - All estimated coefficients for each state.
%                      (AR parameters, standard deviation, transition
%                      matrix - each collum represents each state, ascending order).
%
%   The coefficient std is the model's standard deviation.
%
%   In the Spec_Output.Coeff.p, the row i, collum j, represents the
%   transition probability of being in state j in t-1 and changing to state
%   i on time t.
%
%   In the Coeff.Ar field, the row i, collum j, represents the AR coeff in lag i at
%   state j
%
%   Author: Marcelo Scherer Perlin
%   Email:  [e-mail]
%   Phd Student in finance ICMA/UK 

function [Spec_Output]=MS_AR_Fit(x,ar,k,advOpt)

[nr,nc]=size(x);

if nargin<3
    error('The function needs at least 3 arguments');
end

if nargin==3
    distrib='Normal';
    std_method=1;
    useMex=0;
else
    if isfield(advOpt,'distrib')==0
        distrib='Normal';
    else
        distrib=advOpt.distrib;
    end
    
    if isfield(advOpt,'std_method')==0
        std_method=1;
    else
        std_method=advOpt.std_method;
    end

    if isfield(advOpt,'useMex')==0
        useMex=0;
    else
        useMex=advOpt.useMex;
    end    
    
end

if useMex~=1&&useMex~=0
    error('The input advOpt.useMex only take values 1 or 0');
end

if ~any(std_method==[1 2 3 4])
    error('The input advOpt.std_method should be 1, 2, 3 or 4.');
end

if strcmp(distrib,'Normal')==0&&strcmp(distrib,'t')==0
    error('The distrib input should be ''Normal'' or ''t''');
end

if k==1
    error('If youre using k=1 (one state), you dont need this routine. Just use arma/garch toolbox');
end

if nc>1
    error('The dependent variable should be a vector Tx1.');
end

if ar<1
    error('ar input should be higher than 1')
end

if nr<ar
    error('Your model doesnt make any sense, more lags than observations ??');
end

if nr<=20
    error('You need more than 20 observations in order to estimate the model.')
end

%   Setting up the param0 by estimating a simple Ar(p) model with garch toolbox.
%   If no toolbox is found, the function uses some ordinary guesses.

try
    
    spec=garchset('R',ar,'display','off');
    spec=garchfit(spec,x);
    
    param0_AR=[];
    param0_C=[];
    for i=0:k-1
        param0_AR=[spec.AR -param0_AR];  % Building Ar part for param0 (the next state is always with negative sign from the other)
        param0_C=[spec.C -param0_C];
    end

    switch distrib
        case 'Normal'
            param0=[repmat(std(x),1,k) param0_C param0_AR (nonzeros(repmat(.1,k,k)+eye(k)*(1-k*.1)))'];
        case 't'
            param0=[repmat(std(x),1,k) repmat(5,1,k) param0_C param0_AR (nonzeros(repmat(.1,k,k)+eye(k)*(1-k*.1)))'];
    end

catch

    spec.AR=repmat(.2,1,ar);
    spec.C=mean(x);

    param0_AR=[];
    param0_C=[];
    for i=0:k-1
        param0_AR=[spec.AR -param0_AR];  % Building Ar part for param0 (the next state is always with negative sign from the other)
        param0_C=[spec.C -param0_C];
    end
    
    switch distrib
        case 'Normal'
            param0=[repmat(std(x),1,k) param0_C param0_AR (nonzeros(repmat(.1,k,k)+eye(k)*(1-k*.1)))'];
        case 't'
            param0=[repmat(std(x),1,k) repmat(2,1,k) param0_C param0_AR (nonzeros(repmat(.1,k,k)+eye(k)*(1-k*.1)))'];
    end
end

% Lower and upper bounds at MS estimation

switch distrib
    case 'Normal'
        lB=[repmat(0,1,k)   repmat( -inf,1,k) repmat(-2,1,ar*k) repmat(0,1,k^2)];
        uB=[repmat(inf,1,k) repmat(  inf,1,k) repmat( 2,1,ar*k) repmat(1,1,k^2)];
    case 't'
        lB=[repmat(0,1,k)   repmat(0,1,k)   repmat( -inf,1,k) repmat(-2,1,ar*k) repmat(0,1,k^2)];
        uB=[repmat(inf,1,k) repmat(inf,1,k) repmat(  inf,1,k) repmat( 2,1,ar*k) repmat(1,1,k^2)];
end

warning('off');

global k_global;    % Im using this global here because I need to pass the value of k
% to @confuneq and i haven't found other way around
% it, since such function should have only one input
% parameter

k_global=k;

% Estimation using constrained minimization with @confuneq

options=optimset('fmincon');
options=optimset(options,'display','off');

[param]=fmincon(@(param)MS_AR_Lik(x,param,ar,k,distrib,1,useMex),param0,[],[],[],[],lB,uB,@confuneq_MS_AR,options);

[V]=getvarMatrix_MS_AR(x,param,ar,k,distrib,std_method,useMex);

param_std=sqrt(diag((V))); % the std vector is the square of the diagonal of V (var-cov matrix)

% Controls for covariance matrix. If found imaginary number for variance, replace with
% Inf. This will then be showed at output

param_std(isinf(param_std))=0;

if ~isreal(param_std)
    for i=1:numel(param)
        if ~isreal(param_std(i))
            param_std(i)=Inf;
        end
    end
end

% After finding param, filter it to the data to get estimated output

[arg1,Spec_Output]=MS_AR_Lik(x,param,ar,k,distrib,0,useMex);

% Clearing global variable

clear global k_global;

% calculating the smothed probabilities

Prob_t_1=zeros(nr,k);
Prob_t_1(1,1:k)=1/k; % This is the matrix with probability of s(t)=j conditional on the information in t-1

for i=2:nr
    Prob_t_1(i,1:k)=(Spec_Output.Coeff.p*Spec_Output.filtProb(i-1,1:k)')';
end

filtProb=Spec_Output.filtProb;

P=abs(Spec_Output.Coeff.p);

smoothProb=zeros(nr,k);
smoothProb(nr,1:k)=Spec_Output.filtProb(nr,:);  % last observation for starting filter

smoothProb2=zeros(nr,k);
smoothProb2(nr,1:k)=Spec_Output.filtProb(nr,:);  % last observation for starting filter


for i=nr-1:-1:1     % work backwards in time for smoothed probs
    
    for j1=1:k
        for j2=1:k
            
            smooth_value(1,j2)=smoothProb(i+1,j2)*filtProb(i,j1)*P(j2,j1)/Prob_t_1(i+1,j2);
            
        end
        smoothProb(i,j1)=sum(smooth_value);
    end
    
end

% Calculating Expected Duration of regimes

stateDur=1./(1-diag(Spec_Output.Coeff.p));
Spec_Output.stateDur=stateDur;

% Plotting time varying probabilities

for i=1:k
    States{i}=['State ',num2str(i)];
end

figure(1)
plot([Spec_Output.filtProb]);
xlabel('Time');
ylabel('Filtered States Probabilities');
legend(States);

figure(2)
plot([smoothProb]);
xlabel('Time');
ylabel('Smoothed States Probabilities');
legend(States);

% Sending output to matlab's screen

Spec_Output.smoothProb=smoothProb;
Spec_Output.Coeff.Std_std(1,1:k)=param_std(1:k);

switch distrib
    case 't'
        Spec_Output.Coeff.v_std(1,1:k)=param_std(k+1:k+k);
        Spec_Output.Coeff.const_std(1,1:k)=param_std(k+k+1:k+2*k);
        for i=0:k-1
            Spec_Output.Coeff.AR_std(1:ar,i+1)=param_std(k+2*k+1+i*ar:k+k+k+ar+i*ar);
        end

        for i=0:k-1
            Spec_Output.Coeff.p_std(i+1,1:k)=param_std(end-k^2+1+i*k:end-k^2+k+i*k);
        end
    case 'Normal'
        Spec_Output.Coeff.const_std(1,1:k)=param_std(k+1:2*k);

        for i=0:k-1
            Spec_Output.Coeff.AR_std(1:ar,i+1)=param_std(2*k+1+i*ar:k+k+ar+i*ar);
        end

        for i=0:k-1
            Spec_Output.Coeff.p_std(i+1,1:k)=param_std(end-k^2+1+i*k:end-k^2+k+i*k);
        end

end

fprintf(1,'\n');
fprintf(1,'\n***** Numerical Optimization for MS AR model Converged *****\n\n');
fprintf(1,['Final log Likelihood: ',num2str(Spec_Output.LL),'\n']);
fprintf(1,['Number of parameters: ',num2str(Spec_Output.Number_Parameters),'\n']);
fprintf(1,['Distribution Assumption -> ',distrib,'\n']);
fprintf(1,['Method for standard error calculation -> ',num2str(std_method),'\n']);

fprintf(1,'\n-----> Final Parameters <-----\n');

for j=1:k
    fprintf(1,['\nParameters in State ',num2str(j),' (Coeff value (standard error, p value)): \n\n']);
    fprintf(1,['Residue Standard Deviation: ',num2str(Spec_Output.Coeff.Std(1,j),'%4.4f')]);
    fprintf([' (',num2str(Spec_Output.Coeff.Std_std(1,j),'%4.4f'),', ', ...
                 num2str(2*(1-tcdf(abs(Spec_Output.Coeff.Std(1,j))/Spec_Output.Coeff.Std_std(1,j),nr-numel(param))),'%4.2f'),')\n'])
    fprintf(1,['Mean Constant:              ',num2str(Spec_Output.Coeff.const(1,j),'%4.4f')]);
    fprintf([' (',num2str(Spec_Output.Coeff.const_std(1,j),'%4.4f'),', ', ...
                 num2str(2*(1-tcdf(abs(Spec_Output.Coeff.const(1,j))/Spec_Output.Coeff.const_std(1,j),nr-numel(param))),'%4.2f'),')\n'])

    switch distrib
        case 't'
            fprintf(1,['Degree of Freedom:          ',num2str(Spec_Output.Coeff.v(1,j),'%4.4f')]);
            fprintf([' (',num2str(Spec_Output.Coeff.v_std(1,j),'%4.4f'),', ', ...
                     num2str(2*(1-tcdf(abs(Spec_Output.Coeff.v(1,j))/Spec_Output.Coeff.v_std(1,j),nr-numel(param))),'%4.2f'),')\n'])
    end 

    fprintf(1,'AR parameters:\n');
    for i=1:ar
        fprintf(['     Lag ',num2str(i),': ',num2str(Spec_Output.Coeff.AR(i,j),'%4.4f')]);
        fprintf([' (',num2str(Spec_Output.Coeff.AR_std(i,j),'%4.4f'),', ', ...
                 num2str(2*(1-tcdf(abs(Spec_Output.Coeff.AR(i,j))/Spec_Output.Coeff.AR_std(i,j),nr-numel(param))),'%4.2f'),')\n']);
    end

end

fprintf(1,'\n---> Transition Probabilities Matrix (std. errors in parenthesis) <---\n');

for i1=1:k
    fprintf(1,'\n');
    for i2=1:k
        fprintf(1,['%4.2f (%4.2f)   '],Spec_Output.Coeff.p(i1,i2),Spec_Output.Coeff.p_std(i1,i2));
    end
end

fprintf(1,'\n\n---> Expected Duration of Regimes <---\n\n');

for i=1:k
    fprintf(1,['Expected duration of Regime #%i: %4.2f time periods\n'],i,stateDur(i));
end