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
运行效果如下：
