stvar.m

计量工具箱

for k = 3:nargin
   dk = [ds,num2str(k-2)];
   [row,col] = size(eval(dk));
   if max(row,col) == 1
      scalard(k-2) = 1;
   end
   if row == 1 & scalard(k-2) == 0
      eval([dk,'=',dk,'(:);']);
     row = col;
   end
   if k == 3
      pstring = [lp,dk];
     n = row;
     if nargin == 3
        pstring = [pstring,rp];
      end
   elseif k == nargin & nargin > 3
      pstring = [pstring,c,dk,rp];
   else
      pstring = [pstring,c,dk];
   end
end

% Create index matrix of bootstrap samples.
bootsam = unidrnd(n,n,nboot);

% Get result of bootfun on actual data and find its size. 
thetafit = eval([bootfun,pstring]);
[ntheta ptheta] = size(thetafit);

% Initialize a matrix to contain the results of all the bootstrap calculations.
bootstat = zeros(nboot,ntheta*ptheta);

dbs = 'db';

% Do bootfun - nboot times.
for bootiter = 1:nboot

   for k = 3:nargin
      dk  = [ds,num2str(k-2)];
      dbk = [dbs,num2str(k-2)];
     if scalard(k-2) == 0
         eval([dbk,'=',dk,'(bootsam(:,',num2str(bootiter),'),:);']);
     else
         eval([dbk,'=',dk,';']);
     end
     
      if k == 3
         pstring = [lp,dbk];
        n = row;
        if nargin == 3
           pstring = [pstring,rp];
         end
      elseif k == nargin & nargin > 3
         pstring = [pstring,c,dbk,rp];
      else
         pstring = [pstring,c,dbk];
      end
   end

   evalstr = [bootfun,pstring];
   tmp = eval(evalstr);
   bootstat(bootiter,:) = (tmp(:))';
end

function e = equal(x)

e=x;

function results=stvar_tests(y,param,x)

% PURPOSE: performs a smooth transition vector autoregression
% only presenings tests
%---------------------------------------------------
% USAGE:  result = stvar2(y, param, x)
% where:    y    = an (nobs x neqs) matrix of y-vectors
%        
%           y should be fixed from endog to exog
%
%           Param is a 1x9 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 = the period we want to atart with        
%      transilag = the lag of the transition variable 
%          
%          param = [nlag shockv trans thresh smooth shock sim IRper history transilag]
%
%           (e.g)-> param=[1 2 1 0 1 1 1000 24 30 1]
%
%           x    = optional matrix of variables (nobs x nx)
%                 (NOTE: constant vector automatically included)
%
%       NOTE: Smooth transition function is not applied to exogenous
%             variables
%
%
%---------------------------------------------------
% RETURNS a structure
% results.meth = 'STVAR'
% results.nobs = nobs, # of observations
% results.neqs = neqs, # of equations
% results.nlag = nlag, # of lags
% results.nvar = nlag*neqs+nx+1, # of variables per equation
% --- the following are referenced by equation # --- 
% results(eq).beta  = bhat for equation eq
% results(eq).tstat = t-statistics 
% results(eq).tprob = t-probabilities
% results(eq).resid = residuals 
% results(eq).yhat  = predicted values 
% results(eq).y     = actual values 
% results(eq).sige  = e'e/(n-k)
% results(eq).rsqr  = r-squared
% results(eq).rbar  = r-squared adjusted
% results(eq).boxq  = Box Q-statistics
% results(eq).ftest = Granger F-tests
% results(eq).fprob = Granger marginal probabilities
%---------------------------------------------------
% SEE ALSO: 
%vare, varf, prt_var, prt_granger, prt_ftests (from LeSage's econometrics
% toolbox)
%
%stvar2, stvar_search (from Saki's Toolkit)
%
%---------------------------------------------------

% written by:
% Saki Bigio, Macroeconomic Analysis Dept.
% Banco Central de Reserva del Peru
% Paul de Beaudiez 530,
% Lima L27,  PERU
% [e-mail]

results.meth='STVAR_tests';

[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 = nlag+1;
end;

transilag = param(10);
if transilag == 0
    transilag = 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*1; % # of variables per equation

%Building the Smooth Transition Function
transi   = y(:,trans);
standard = std(transi);
transi   = lag(transi,transilag);
gfunc    =(1+exp(-smooth*(transi-thresh)/standard)).^(-1); %Weise uses the function with -1/2
g        = transi;

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

%Building the non-linear data

gy=[]; 
ylag=mlag(y,nlag);

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

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

%We should fix the first period cause of the lag and SMTV
varres=vare2(y,nlag,gy);

%Loading Variables for Impulse Response
n = varres(1).nobs;
k = varres(1).neqs;

%get coefficients %%%%%%Generalize
b=[];
for i=1:neqs          
    b  = [b varres(i).beta];
end

%Extract reduced form residuals and (using Cholesky decomposition, obtain
%functional form uncorrelated errores)
e = zeros(n,k);
for i=1:k
    e(:,i) = varres(i).resid;
end

%compute tests statistics for non-linearity
%start with the equation by equation test
varres_aux=vare2(y,nlag);

res_r=zeros(n,k);
for i=1:k
    res_r(:,i) = varres_aux(i).resid;   %get restricted residuals
end;

res_ur=e;                                %get unrestricted residuls

%Calculating the Equation by Equation LM estimator
FStat_eq=zeros(neqs,1);
for i=1:neqs
    SSR0=sum(res_r(:,i).^2);
    SSR1=sum(res_ur(:,i).^2);
    FStat_eq(i)=((SSR0-SSR1)/(nlag*neqs))/((SSR1/((nobs-nlag)-(2*nlag*neqs+1))));%calculate small sample F-Statistic
end;

pval=1-fdis_prb(FStat_eq,nlag*neqs,(nobs-nlag)-(2*nlag*neqs+1));        %calculate corresponding p-value
results.fstat_pval=pval;

%the overall LR-test

%omega0= (1/((nobs-nlag)-neqs*nlag-1))*res_r'*res_r;                  %build cov. matrices
%omega1= (1/((nobs-nlag)-2*neqs*nlag-1))*res_ur'*res_ur;

%LR=((nobs-nlag)-nlag*neqs-1)*(log(det(omega0))-log(det(omega1)));
%LRpval=1-chis_prb(LR,nlag*neqs^2);

omega0= (1/((nobs-nlag)))*res_r'*res_r;                  %build cov. matrices
omega1= (1/((nobs-nlag)))*res_ur'*res_ur;

LR=(nobs-nlag)*(log(det(omega0))-log(det(omega1)));
LRpval=1-chis_prb(LR,(nlag*neqs^2));

results.omega0=omega0;
results.omega1=omega1;      
results.LR=LR;
results.LRpval=LRpval;
results.beta=b;
results.g=g;


function results = vare2(y,nlag,x)
% PURPOSE: performs vector autogressive estimation and presents no tests
%---------------------------------------------------
% USAGE:  result = vare(y,nlag,x) 
% where:    y    = an (nobs x neqs) matrix of y-vectors
%           nlag = the lag length
%           x    = optional matrix of variables (nobs x nx)
%                 (NOTE: constant vector automatically included)
%---------------------------------------------------
% RETURNS a structure
% results.meth = 'vare'
% results.nobs = nobs, # of observations
% results.neqs = neqs, # of equations
% results.nlag = nlag, # of lags
% results.nvar = nlag*neqs+nx+1, # of variables per equation
% --- the following are referenced by equation # --- 
% results(eq).beta  = bhat for equation eq
% results(eq).tstat = t-statistics 
% results(eq).tprob = t-probabilities
% results(eq).resid = residuals 
% results(eq).yhat  = predicted values 
% results(eq).y     = actual values 
% results(eq).sige  = e'e/(n-k)
% results(eq).rsqr  = r-squared
% results(eq).rbar  = r-squared adjusted
% results(eq).boxq  = Box Q-statistics
% results(eq).ftest = Granger F-tests
% results(eq).fprob = Granger marginal probabilities
%---------------------------------------------------
% SEE ALSO: varf, prt_var, prt_granger, prt_ftests 
%---------------------------------------------------

% written by:
% James P. LeSage, Dept of Economics
% University of Toledo
% 2801 W. Bancroft St,
% Toledo, OH 43606
% [e-mail]

[nobs neqs] = size(y);

results(1).meth = 'vare2';

nx = 0;

if nargin == 3
[nobs2 nx] = size(x);
 if (nobs2 ~= nobs)
 error('vare2: nobs in x-matrix not the same as y-matrix');
 end;
end;

% adjust nobs to feed the lags
nobse = nobs - nlag;

% nvar adjusted for constant term 
 k = neqs*nlag + 1 + nx;
 nvar = k;

results(1).nvar = nvar;

xlag = mlag(y,nlag);

results(1).nobs = nobse;
results(1).neqs = neqs;
results(1).nlag = nlag;


% form x-matrix
if nx 
xmat = [xlag(nlag+1:nobs,:) x(nlag+1:nobs,:) ones(nobs-nlag,1)];
else
xmat = [xlag(nlag+1:nobs,:) ones(nobs-nlag,1)];
end;


% pull out each y-vector and run regressions
for j=1:neqs;

 yvec = y(nlag+1:nobs,j);
 res = ols(yvec,xmat);
 results(j).beta  = res.beta;      % bhats
 results(j).tstat = res.tstat;     % t-stats
 % compute t-probs
      tstat = zeros(nvar,1);
      tstat = res.tstat;
      tout = tdis_prb(tstat,nobse-nvar);
 results(j).tprob = tout;          % t-probs
 results(j).resid = res.resid;     % resids 
    sigu = res.resid'*res.resid;
 results(j).yhat = res.yhat;       % yhats
   results(j).y    = yvec;           % actual y
   results(j).rsqr = res.rsqr;       % r-squared
   results(j).rbar = res.rbar;       % r-adjusted
   results(j).sige = res.sige;


% do the Q-statistics
% use residuals to do Box-Pierce Q-stats
% use lags = nlag in the VAR
% NOTE: a rule of thumb is to use (1/6)*nobs
%       but this seems excessive to me
elag = mlag(res.resid,nlag);
% feed the lags
etrunc = elag(nlag+1:nobse,:);
rtrunc = res.resid(nlag+1:nobse,1);
qres   = ols(rtrunc,etrunc);
if nlag ~= 1
 boxq   = (qres.rsqr/(nlag-1))/((1-qres.rsqr)/(nobse-nlag));
else
 boxq   = (qres.rsqr/(nlag))/((1-qres.rsqr)/(nobse-nlag));
end;

results(j).boxq = boxq;

% form x matrices for joint F-tests
% exclude each variable from the model sequentially

for r=1:neqs;
xtmp = [];
for s=1:neqs;
 if s ~= r
   xlag = mlag(y(:,s),nlag);
 xtmp = [xtmp trimr(xlag,nlag,0)];
   end;
end;
% we have an xtmp matrix that excludes 1 variable
% add deterministic variables (if any) and constant term
if nx > 0
xtmp = [xtmp x(1:nobse,:) ones(nobse,1)];
else
xtmp = [xtmp ones(nobse,1)];
end;
% get ols residual vector
b = xtmp\yvec; % using Cholesky solution
etmp = yvec-xtmp*b;
sigr = etmp'*etmp;
% joint F-test for variables r
ftest(r,1) = ((sigr - sigu)/nlag)/(sigu/(nobse-k)); 
end;

results(j).ftest = ftest;     
%results(j).fprob = fdis_prb(ftest,nlag,nobse-k);

end; 
% end of loop over equations