stvar.m

计量工具箱

function results = stvar(y, param, x)
% PURPOSE: performs a smooth transition vector autoregression
%---------------------------------------------------
% USAGE:  result = stvar(y, param, x)
% where:    y    = an (nobs x neqs) matrix of y-vectors
%        
%           y should be fixed from most endogenous to most exogenous
%
%           Param is a 1x10 vector that includes the following 
%           information in order:
%
%           nlag = the lag length
%         shockv = variable of y being shocked (column of y)
%          trans = position of the transition variable in the y matrix
%         thresh = the threshold value
%         smooth = the smoothnes parameter                
%          shock = the standard deviation shocks 
%            sim = the number of montecarlo simulations
%          IRper = the number of periods for the impulse response function
%        history = 0 if choosen period under threshold (1 if not)
%       translag = the lag of the transition variable  
%
%          param = [nlag shockv trans thresh smooth shock sim IRper history translag]
%
%           (e.g.)-> param=[1 2 1 0 1 1 1000 24 1 1]
%
%           x    = optional matrix of variables (nobs x nx)
%                 (NOTE: constant vector automatically included)
%
%       NOTE: Smooth transition function is not applied to exogenous
%             variables
%
%---------------------------------------------------
%
%      Bibliography:
%
%      The non linear VAR is built on the following structure (reduced form): 
%
%      Y(t)=[(A-B(g(t))L]Y(t)+E(t) 
%
%      where  
%
%      g(t)=is the logistic form depending on a state variable
%
%      and L is a Lag Operator
%
%       For furtherings check out Wiese 1999
%---------------------------------------------------
% RETURNS a structure
% results.meth = 'STVAR'
% results.nobs = nobs, # of observations
% results.neqs = neqs, # of equations
% results.nlag = nlag, # of lags
%
% Generalized Impulse Response Function (see Pesaran, Potter & Koop (1997))
%
% results.irfs      = Generalized Impulse Response Function (GIRF)
% results.irfs_sup  = Superior Band of the GIRF
% results.irfs_inf  = Inferior Band of the GIRF
%
% Non-Linearity Tests (see Terasvirta (1995))
%
% results.omega0    = Var-Cov of the equivalent linear VAR
% results.omega1    = Var-Cov;
% results.LR        = tests.LR;
% results.LRpval    = tests.LRpval;
%
% take into account that results (for parameters and coefficients are presented in the following order):   
% 
% y1(t-1), y1(t-2),....y1(t-nlag),y2(t-1)....ynvar(t-nlag),y1(t-1),
% g*y1(t-2),....g*y1(t-nlag),g*y2(t-1)....g*ynvar(t-nlag)
%
%---------------------------------------------------
% SEE ALSO: vare, varf, prt_var, prt_granger, prt_ftests (from LeSage's econometrics
% toolbox) and optsmooth (from Saki's toolkit)
%
%---------------------------------------------------
%
% written by:
% Saki Bigio 
% Department of Macroeconomic Analysis, 
% Banco Central de Reserva del Peru
% Paul de Beaudiez 530,
% Lima L27,  PERU
% [e-mail]

results.meth='Smooth Transition VAR';

[nobs neqs] = size(y);

results.nobs = nobs; % # of observations
results.neqs = neqs; % # of equations

%Setting Default Values
nlag   = param(1);
if nlag ==0
    nlag=1;
end;

results.nlag = nlag; % # of lags

shockv   = param(2);
if shockv == 0
    shockv = neqs;
end;

trans    = param(3);
if trans == 0
    trans = 1;
end;

thresh   = param(4);
if thresh == 0
    thresh = 0;
end;

smooth   = param(5);
if smooth == 0
    smooth = 7;
end;

shock    = param(6);
if shock == 0
    shock = 1;
end;

sim      = param(7);
if sim == 0
    sim   = 1000;
end;

IRper    = param(8);
if IRper == 0
    IRper = 24;
end;

Hist     = param(9);
if Hist == 0 | Hist == 1
else
    error('choose a type of history');
end;

translag = param(10);
if translag == 0
    translag = 1;
end;

if nargin == 3
 [nobsx nx] = size(x);
 if (nobsx ~= nobs)
  error('var: nobs in x-matrix not the same as y-matrix');
 end;
results.nvar = 2*nlag*neqs+nx+2; % # of variables per equation
end;

results.nvar = 2*nlag*neqs+2;    % # of variables per equation

%Building the Smooth Transition Function
transi = y(:,trans);
standard=std(transi);
transi = lag(transi,translag);
gfunc = (1+exp(-smooth*(transi-thresh)/standard)).^(-1); %Weise uses the function with -1/2, note that variables aren't normalized
g=gfunc;

%Generalizing nx for case when we have exog. var.
nx = 0;

%Building the non-linear data matrix
gy=[]; 
ylag=mlag(y,nlag);

for i=1:neqs*nlag;
gy=[gy g.*ylag(:,i) ];
end    

if nargin == 3
gy=[gy x];
end;  

%Running the Non-Linear VAR 
varres=vare2(y,nlag,gy);

%Loading Variables for Impulse Response
n = varres(1).nobs;
k = varres(1).neqs;
gy=gy((1+nlag:nobs),:);
g=g((1+nlag:nobs),:);

%Estimation Coefficients 
b=[];
for i=1:neqs          
    b  = [b varres(i).beta];
end

%Estimation T-Stats
tstat=[];
for i=1:neqs          
    tstat  = [tstat varres(i).tstat];
end

%Estimation P-Values
tprob=[];
for i=1:neqs          
    tprob  = [tprob varres(i).tprob];
end

%Estimation Residuals
e = zeros(n,k);
for i=1:k
    e(:,i) = varres(i).resid;
end

% Recovering Predicted Values
yhat=[];
for i=1:neqs          
    yhat  = [yhat varres(i).yhat];
end

% Recovering Alligned Actual Values
yact=[];
for i=1:neqs          
    yact  = [yact varres(i).y];
end

% Recovering R-squared
rsqr=[];
for i=1:neqs          
    rsqr  = [rsqr varres(i).rsqr];
end

% Recovering Alligned Actual Values
yact=[];
for i=1:neqs          
    yact  = [yact varres(i).y];
end

% Recovering R-Squared Adjusted
rbar=[];
for i=1:neqs          
    rbar  = [rbar varres(i).rbar];
end

%registering relevant statistics of the estimation
results.beta  = b;
results.tstat = tstat;
results.tprob = tprob;
results.resid = e;
results.yhat  = yhat ;
results.y     = yact;
results.rsqr  = rsqr;
results.rbar  = rbar;

%Extract reduced form residuals and (using Cholesky decomposition, obtain
%functional form uncorrelated errores)
Omega = (1/(n-k))*e'*e;
C = chol(Omega);
v = inv(C)*e';
v=v';

%Creating empty matrix to store differenced paths
paths_dif=[];

%Preparing x
if nargin == 3
        x=x((nlag+1:nobs),:);    
end

%Starting main simulation loop
for i = 1:sim

    %Choosing a Random period that coincides with the relevant history
    junk=randint(1,1,nlag+1, nobs-2*IRper);
    
    if Hist ==1
        while y(junk-translag,trans) < thresh
              junk=randint(1,1,nlag+1, nobs-2*IRper);
        end
        History=junk;
    else
        while y(junk-translag,trans) > thresh
              junk=randint(1,1,nlag+1, nobs-2*IRper);
        end
        History=junk;
    end
        
    %Reshufling uncorrelated errors sim times (montecarlo simulations)
    %and fixing size

    temp = bootstrp(1,'equal',v);
    vbs = reshape(temp,n,k);
    
    %Fix a given shock because we are using Cholesky decomposition, a direct addition can be
    %done that will be equal to fixing a "shock" sized standard deviation shock
    %(see Hamilton p. 400)
    
    vbss=vbs;
    vbss(History,shockv)= shock; %???this is only the shock, should we add it to the residual...?  
    vbs(History,shockv)= 0;
    
    %returning dependence to shocks

    ebs_  = (C*vbs')';
    ebss_ = (C*vbss')';
    eb    = [ebs_ ebss_];

    %Setting an auxiliary variable and all complements to rebuild series
    paths=[];
    for i=0:neqs:neqs
    
    
        ebs=eb(:,i+1:i+neqs);
        path = y;
        g_sim=gfunc((1+nlag:nobs),:);
    
        ylag = mlag(y,nlag);
        ymat=ylag(nlag+1:nobs,:);

        gy=[];
        for i=1:neqs*nlag   
            gy=[gy g_sim.* ymat(:,i)]; 
        end   
    
        if nargin == 3
            gy=[gy x];    
        end; 
       
        gy=[ymat gy ones(nobs-nlag,1)];

        %Rebuilding series for I-R period
        for t=History:History+IRper; % (le quitamos 
            path(t,:)= gy(t,:)*b + ebs(t,:);

            %updating gy matrix with the values just obtained
            for j = 0:k-1
                ymat(t+1,j*nlag+1:(j+1)*nlag)=rot90(path((t-nlag+1):t,j+1),-1); %check this out (t or t+1)
            end;

            %Updating point t of the G function 
            g_sim(t) = (1+exp(-smooth*(path(t-translag,trans)-thresh)/standard))^(-1);

            %Emptying GY matrix and updating it
            gy=[];

            for i=1:neqs*nlag
                gy=[gy g_sim.*ymat(:,i)];
            end;    
        
            if nargin == 3
                gy=[gy x];    
            end;
    
            gy=[ymat gy ones(nobs-nlag,1)]; 
    
        end;
    
        path=path((History):(History+IRper),:); 
        %store resulting path
        paths=[paths path];
    
    end;

    paths_us  = paths(:,1:neqs);
    paths_s   = paths(:,neqs+1:2*neqs);
    paths_dif = [paths_dif paths_s-paths_us];

end;  

%Reordering main loop's output

for j=1:neqs
    for i=1:sim
    IRFS_aux(:,i+(j-1)*sim)=paths_dif(:,j+neqs*(i-1));
    end;
end;


%Building I-R and confidence bands over the 90%
inf(IRper+1,neqs)=0;
sup(IRper+1,neqs)=0;

for i=1:neqs
    for t=1:IRper+1
        IRF(t,i)=mean(IRFS_aux(t,(1+(i-1)*sim):i*sim));
        [prob val]=ecdf(IRFS_aux(t,(1+(i-1)*sim):i*sim));
        inf(t,i)=val(min(find(prob>0.05)));
        sup(t,i)=val(min(find(prob>=0.95)-1));
    end;
end;

results.irfs=IRF;
results.irfs_sup=sup;
results.irfs_inf=inf;

%Graphing Impulse Response Function
periods=1:IRper+1;

for j=1:neqs
    subplot(neqs,1,j);
    plot(periods,sup(:,j),'r',periods,inf(:,j),'r',periods,IRF(:,j),'b-',periods,0,'k');
    xlabel('Periods');
    ylabel('variable ?');
end

%compute tests statistics for non-linearity using stvar_tests
tests=stvar_tests(y,param);

results.fstat_pval=tests.fstat_pval;
results.omega0=tests.omega0;
results.omega1=tests.omega1;
results.LR=tests.LR;
results.LRpval=tests.LRpval;

function a = randint(m,n,a,b)
% RANDINT Randomly generated integral matrix
%	randint(m,n) returns an m-by-n matrix with entries between
%	0 and 9.
%	randint(m,n,a,b) returns entries between a and b.
if nargin < 3, a = 0; b = 9; end
a = floor((b-a+1)*rand(m,n))+ a;

function [bootstat, bootsam] = bootstrp(nboot,bootfun,d1,d2,d3,d4,d5,d6,d7,d8,d9,d10);
%BOOTSTRP Bootstrap statistics.
%   BOOTSTRP(NBOOT,BOOTFUN,D1,...) draws NBOOT bootstrap data samples and
%   analyzes them using the function, BOOTFUN. NBOOT must be a positive integer.
%   BOOTSTRAP passes the (data) D1, D2, etc. to BOOTFUN.
%
%   [BOOTSTAT,BOOTSAM] = BOOTSTRP(...) Each row of BOOTSTAT contains
%   the results of BOOTFUN on one bootstrap sample. If BOOTFUN returns a matrix,
%   then this output is converted to a long vector for storage in BOOTSTAT.
%   BOOTSAM is a matrix of indices into the row

%   Reference:
%      Efron, Bradley, & Tibshirani, Robert, J.
%      "An Introduction to the Bootstrap", 
%      Chapman and Hall, New York. 1993.

%   B.A. Jones 9-27-95
%   Copyright (c) 1993-98 by The MathWorks, Inc.
%   $Revision: 2.5 $  $Date: 1997/11/29 01:44:59 $

% Initialize matrix to identify scalar arguments to bootfun.
scalard = zeros(nargin-2,1);

lp = '(';      % Left parenthesis string variable.
rp = ')';      % Right parenthesis string variable.
c  = ',';      % Comma string variable.
ds = 'd';      % 'd' as a string variable.

% Construct argument list for bootfun