VAR模型的matlab代码

var模型的matlab代码，摘自：https://sites.google.com/site/jnakajimaweb/var

%%--------------------------------------------------------%%
%%           Estimation of standard VAR models            %%
%%--------------------------------------------------------%%
%%              coded by Jouchi Nakajima                  %%
%%--------------------------------------------------------%%
%%               Last update: 2018/07/31                  %%
%%       http://sites.google.com/site/jnakajimaweb/       %%
%%--------------------------------------------------------%%
%%    You may use and modify this code at your own risk   %%
%%--------------------------------------------------------%%

%%  OLS estimation of standard (time-invariant parameter)
%%  VAR models with Cholesky-type identification to
%%  produce impulse response and historical decomposition

clear all
close all

%% ------------------- load data ------------------------ %%

% specify file and sheet
mdata = xlsread('test.xlsx', 'SampleData'); 

% specify variable names
asvar = {'Inf', 'X', 'r', 'Q'};

%% ------------------- set options ---------------------- %%

%%%% model and series settings

vind = [1 2 3 4];   % select column index of data
% vind = [4 2 1 3]; % (eg: different ordering)
% vind = [3 2 1];   % (eg: sub-sets of data)

flc = 1;         % constant term (1:yes, 0:no)

nl = 1;          % # of lags

%%%% impulse response

vsign = [1 1 1 1];     % signs of shocks
% vsign = [1 -1 -1 1]; % (eg: negative shocks for some variables)

flsu = 0;        % shock size (1:one unit, 0:one standard error)

dband = 1.96;    % width of condidence intervals
                 % (one standard error times "dbands")

vac = [0 0 0 0];   % accumulated response (1: yes, 0: no)
% vac = [0 1 0 1]; % (eg: accumulated for some variables)

nimp = 60;       % # of horizons to draw

flcum = 1;       % label for accumulated responses (1:on, 0:off)

%%%% historical decomposition

vhd = [0 1 1 0]; % historical decomposition (1:yes, 0:no)

%%%% output files for figures

flout = 1;       % output file (1: yes, 0: no)

%%%% simulation settings

nsim = 10^3;     % simulation size for confidence intervals
                 % (0: no intervals)

irs = 7;         % random seed

%% ------------------------------------------------------ %%

%% set variables %%

my = mdata(:, vind);    % select data
asvar = asvar(vind);

[ns, nk] = size(my);    % # of times, variables
nsl = ns - nl;          % # of times to estimate
nb = flc + nk * nl;     % # of coefficients

rand('seed', irs);      % random seed
randn('seed', irs);

if size(vsign, 2) ~= nk
  fprintf('\n%%%% ! The size of "vsign" should equal')
  fprintf('\n%%%%   the number of variables\n')
end
if size(vac, 2) ~= nk
  fprintf('\n%%%% ! The size of "vac" should equal')
  fprintf('\n%%%%   the number of variables\n')
end
if sum(vhd) > 0 && size(vhd, 2) ~= nk
  fprintf('\n%%%% ! The size of "vhd" should equal')
  fprintf('\n%%%%   the number of variables\n')
end
if size(vsign, 2) > nk
  vsign = vsign(1:nk);
end
if size(vac, 2) > nk
  vac = vac(1:nk);
end
if sum(vhd) > 0 && size(vhd, 2) > nk
  vhd = vhd(1:nk);
end


%% estimate VAR %%

fprintf('\n%%%% [Estimate VAR]')

mx = [ones(ns-nl, flc), zeros(ns-nl, nk*nl)];
for j = 1 : nl
    mx(:, (j-1)*nk+flc+1:j*nk+flc) = my(nl-j+1:ns-j,:);
end
mx2i = inv(mx' * mx);

mb = zeros(nb, nk);
me = zeros(nsl, nk);
vl = zeros(nk, 1);

for i = 1 : nk
   vy = my(nl+1:end, i);

   vb = mx2i * (mx' * vy);
   ve = vy - mx * vb;
   de2 = ve' * ve / nsl;

   vl(i) = -0.5 * nsl * (log(2 * pi * de2) + 1);

   mb(:, i) = vb;
   me(:, i) = ve;
end

mes = me' * me;
mS = mes / (nsl - nb);
met = (chol(mS, 'lower') \ me')';

vaic = -2 * vl + 2 * nb;
vbic = -2 * vl + 2 * nb * log(nsl);

fprintf('\n%%%%')
fprintf('\n%%%% Lag:%2.0f', nl)
fprintf('\n%%%% Const: ')
if flc == 1, fprintf('on')
else, fprintf('off'), end
fprintf('\n%%%%')
fprintf('\n%%%% AIC: %7.1f', sum(vaic))
fprintf('\n%%%% BIC: %7.1f', sum(vbic))
fprintf('\n%%%%\n')


%% compute impulse response %%

mimp = zeros(nimp+1, nk^2);
mimpm = mimp;
mimps = mimp;

vbh = mb(:);
mW = inv(mes);

mid = reshape((1 : nk*nl), nl, nk)';
vid = mid(:)';

flac = sum(vac) > 0;
vack = kron(ones(1, nk), vac) > 0;

mbi = mb;
mSi = mS;
for ii = 0 : nsim

  if ii > 0
    mSi = inv(wishrnd(mW, nsl+nk));
    mVi = kron(mSi, mx2i);
    mVi = (mVi + mVi') / 2;
    mbi = reshape(mvnrnd(vbh, mVi), nb, nk);
  end
  mbj = mbi(1+flc:end, :);
  vbi = mbj(:);
  mLi = chol(mSi, 'lower');

  for j = 1 : nk

    vshock = zeros(nk, 1);
    vshock(j) = vsign(j);

    for t = 1 : nimp+1
      if t == 1
        myt = zeros(nl, nk);
        myt(1, :) = (mLi * vshock)';
      else
        vyt = myt(:)';
        mxt = kron(eye(nk), vyt(vid));
        vyt1 = (mxt * vbi)';
        if nl > 1
            myt = [vyt1; myt(1:end-1, :)];
        else
            myt = vyt1;
        end
      end
      mimp(t, (j-1)*nk+(1:nk)) = myt(1, :);
    end
  end

  if flac
    mimp(:, vack) = cumsum(mimp(:, vack));
  end

  if ii == 0
    mimph = mimp;
  else
    mimpm = mimpm + mimp;
    mimps = mimps + mimp.^2;
  end

end

if nsim > 0
  mimpm = mimpm / nsim;
  mimps = real(sqrt(mimps / nsim - mimpm.^2));

  mimpu = mimph + mimps * dband;
  mimpl = mimph - mimps * dband;
end

%% plot

if flsu == 1
  ves = sqrt(diag(mS));
else
  ves = ones(nk, 1);
end
if flout == 1
  mout = zeros(nimp+1, nk^2*3);
end
figure
for i = 1 : nk
  for j = 1 : nk
    kk = (i-1)*nk + j;
    subplot(nk, nk, (j-1)*nk+i)
    hold on
    if nsim > 0
      plot(0:nimp, mimph(:, kk) / ves(i), 'b-')
      plot(0:nimp, mimpu(:, kk) / ves(i), 'r--')
      plot(0:nimp, mimpl(:, kk) / ves(i), 'r--')
    else
      plot(mimph(:, kk) / ves(i), '-')
    end
    vax = axis;
    axis([0 nimp vax(3:4)])
    if vax(3) * vax(4) < 0
      line([0, nimp], [0, 0], 'Color', ones(1,3)*0.6)
    end
    if vac(j) == 1 && flcum == 1
      sac = '(cum)';
    else
      sac = '';
    end
    title(['Response of ', char(asvar(j)), sac, ...
           ' to ', char(asvar(i))], 'Fontsize', 9)
    hold off

    if nsim > 0 && flout == 1
      mout(:, (kk-1)*3 + (1:3)) ...
      = [mimpl(:, kk) mimph(:, kk) mimpu(:, kk)] / ves(i);
    end
  end
end

%% output file

if flout == 1
  flb = nsim ~= 0;
  if ~flb
    mout = mimph;
  end
  if isequal(exist('var_imp.xlsx', 'file'), 2)
    delete('var_imp.xlsx');
  end
  
  asl = cell(3, nk^2*3+1);
  asl{1, 1} = 'Response:';
  asl{2, 1} = 'Shock:';
  asb = {'low', 'mean', 'high'};
  ii = 2;
  for i = 1 : nk
  for j = 1 : nk
  for k = 1 : 1 + flb * 2 
    if vac(j) == 1 && flcum == 1
      sac = ' (cum)';
    else
      sac = '';
    end
    if k == 1
      asl{1, ii} = [asvar{j}, sac];
      asl{2, ii} = asvar{i};
    end
    asl{3, ii} = asb{k+~flb};
    ii = ii + 1;
  end
  end
  end
  
  xlswrite('var_imp.xlsx', asl, 'Sheet1', 'A1');
  xlswrite('var_imp.xlsx', [(0:nimp)', mout], 'Sheet1', 'A4');
end


%% compute historical decomposition %%

if sum(vhd) > 0

 if flc == 1
   vc = mb(1, :)';
 end

 mbj = mb(1+flc:end, :);
 vbi = mbj(:);
 mLi = chol(mS, 'lower');

 myh = zeros(ns, (nk+flc)*nk);

 for i = 1-flc : nk
  myt = zeros(nl, nk);
  for t = nl+1 : ns
    if ~i
      vet = vc;
    else
      vshock = zeros(nk, 1);
      vshock(i) = met(t-nl, i);
      vet = mLi * vshock;
    end

    vyt = myt(:)';
    mxt = kron(eye(nk), vyt(vid));

    vyt1 = (mxt * vbi + vet)';
    if nl > 1
      myt = [vyt1; myt(1:end-1, :)];
    else
      myt = vyt1;
    end

    myh(t, (0:nk-1)*(flc+nk)+i+flc) = vyt1;

  end
 end

 asleg = {'Data'};
 if flc == 1
   asleg = [asleg, {'const'}];
 end
 asleg = [asleg, asvar];

 for i = 1 : nk
  if vhd(i) == 1
   figure
   plot(my(:, i), '-')
   hold on
   plot(nl+1:ns, myh(nl+1:end, (1:flc+nk)+(i-1)*(flc+nk)), ...
        '--')
   hold off
   vax = axis;
   axis([1 ns vax(3:4)])
   if vax(3) * vax(4) < 0
    line([1, ns], [0, 0], 'Color', ones(1,3)*0.6)
   end
   legend(asleg)
   title(['Historical decomposition of' asvar(i)], ...
         'FontSize', 12)
  end
 end

 if flout == 1
  if isequal(exist('var_hd.xlsx', 'file'), 2)
    delete('var_hd.xlsx');
  end
  asl = cell(2, nk*(nk+flc)+1);
  asl{1, 1} = 'Explained:';
  asl{2, 1} = 'Shock:';
  ii = 2;
  for i = 1 : nk
  for j = 1 : nk+flc
    asl{1, ii} = asleg{i+2};
    asl{2, ii} = asleg{j+1};
    ii = ii + 1;
  end
  end
  xlswrite('var_hd.xlsx', asl, 'Sheet1', 'A1');
  xlswrite('var_hd.xlsx', [(1:ns)', myh], 'Sheet1', 'A3');
 end

end

运行效果如下：