matlab 后缀名为model,MATLAB SVAR model toolboxs

%% SVAR: An application

% Paper: 'Identifying the interdependence between US monetary policy and

% the stock market', Bjornland and Leitemo (2009)

%

% *************************************************************************

% By Alex Carrasco, november 2017

% *************************************************************************

tic

clear all; clc; close all;

boot = 1; % bootstrap again?

B = 50; % numbers of replications [2000]

K = 100; % total periods for IRF

% set up for printing figures

% path_g = [pwd '\Graphs'];

% imprime = @(x) print( gcf, '-depsc2', [path_g filesep x]);

% imprpdf = @(x) eps2pdf( [path_g filesep x '.eps']);

formataxis = @(x) set(x, 'Fontname', 'Times', 'FontWeight', 'normal', 'Fontsize', 23, 'Box', 'Off', 'PlotBoxAspectRatio', [1 0.75 1]);

formatlegend = @(x) set(legend, 'Location', x, 'Orientation', 'Vertical', 'Box', 'Off', 'Fontsize', 22, 'Fontangle', 'normal');

%% [I] loading dataset and data transformation

% [a] loading

load data;

ran = se(1):47:se(end); % range for plots [do not change!]

m=size(DD,2);

% [b] unit root tests

UR=uroottest(DD,'report','ur_test','names',var_names,'alpha',0.05); % creates a table with test results

% [c] plotting data

cc=1;

for ii=1:m

figure(cc);

h=plot(se,DD(:,ii),'-','linewidth',2.5); axis tight;

set(h,'color',[0 0 0.5]);

x=gca;

set(x,'xtick',ran,'xticklabel',ranlab);

formataxis(x);

% imprime([cod_names{ii} '_data']);

% imprpdf([cod_names{ii} '_data']);

cc=cc+1;

end

%% [II] Model selection, p: FPE, AIC, BIC, HQ

pmax=8;

[popt,ic]=VAROptLag(DD,pmax);

fprintf('\n Optimal lags: %1.0f %1.0f %1.0f %1.0f\n', popt);

%% [III] Reduced form estimation: VAR(p)

p = popt(3); % BIC optimal lag

% [a] Baseline model

[Chat,Shat,F,Theta]=VARest(DD,p);

% [b] Alternative model

DD(:,[4 5])=DD(:,[5 4]);

[Chata,Shata,Fa]=VARest(DD,p);

DD(:,[4 5])=DD(:,[5 4]);

P = chol(Shat); % Cholesky decomposition (baseline)

Pa = chol(Shata); % Cholesky decomposition (alternative)

%% [IV] A-model, identification and estimation

% RA is the restriction matrix

A = nan(m,m);

% [a] Zero restrictions

A(1,2:5) = 0; A(2,3:5)=0; A(3,4:5)=0;

free_param = isnan(A);

nfree = sum(vec(free_param));

ind0 = reshape(1:m^2,m,m);

ind0 = ind0(free_param);

RA = zeros(m^2,nfree);

[I,J] = ind2sub([m,m],ind0);

for k=1:nfree

RA((J(k)-1)*m+I(k),k)=1;

end

% [b] Long-run restriction

RA(m^2,m*(m+1)/2)=-Theta(4,4)/Theta(4,5);

RA(:,m*(m+1)/2+1)=[];

% [c] Estimation:

% Based on Amisano and Gianini (1997)/ Marquardt algorithm (Baseline ordenation)

gamma0=.8*ones(nfree-1,1); gamma0(nfree-1)=-gamma0(nfree-1); % Sign restriction

[gammaA,logL,error,S]=mlSVAR(RA,gamma0,Shat,T,'tau',0.3,'errotol_score',1e-7,'errortol_grad',1e-7);

A=reshape(RA*gammaA,m,m);

%% [V] IRF [mean]

% [a] Structural identifcation

IR=VARimpulse(F,A,K); % mean responses in structural identifcation

% [b] Models identificated with recursive restrictions

Kmax = 46; % toplot

x =(1:Kmax)';

sh = [5 4]; % shocks: i/Ds \\ Ds/i

resp = [5 4]; % shocks: i/Ds \\ Ds/i

IRchol1 =VARimpulse(F,P',K); % IRF using baseline order

IRchol2 =VARimpulse(Fa,Pa',K); % IRF using alternative order

for ss=1:numel(sh)

for vv=1:numel(resp)

figure(cc);

irf1 = IRchol1(resp(vv),:,sh(ss))'/IRchol1(sh(ss),1,sh(ss)); % normalizing responses

irf2 = IRchol2(resp(resp(vv)~=resp),:,sh(sh(ss)~=sh))'/IRchol2(sh(sh(ss)~=sh),1,sh(sh(ss)~=sh));

if resp(vv)==4

irf1=cumsum(irf1);

irf2=cumsum(irf2);

end

irf1=irf1(x); irf2=irf2(x);

h=plot([irf1,irf2]); hold on;

plot([1 Kmax],[0 0],':k');

set(h(1),'color',[0 0 .5],'linewidth',3);

set(h(2),'color',[0 0 .5],'linestyle','--','linewidth',2);

% axis

axs=gca; miny=min(min([irf1,irf2])); maxy=max(max([irf1,irf2]));

xlim([1 Kmax]);

ylim([min(miny,0)*1.2 maxy*1.2]);

set(axs,'xtick',6:5:Kmax,'Xticklabel',5:5:Kmax-1);

xlabel('Months','fontname','Times','FontWeight', 'normal', 'Fontsize', 22);

ylabel('Percent','fontname','Times','FontWeight', 'normal', 'Fontsize', 22);

formataxis(axs);

if resp(vv)==5

lg=legend('Baseline','Alterntive');

formatlegend('NorthEast');

end

% imprime(['IRFChol_' cod_names{resp(vv)} '_' cod_names{sh(ss)}]);

% imprpdf(['IRFChol_' cod_names{resp(vv)} '_' cod_names{sh(ss)}]);

cc=cc+1;

end

end

%% [VI] Standard errors [Bootstraping]

if boot

[bhat,Omg]=VAR_boots(DD,'order',p,'n_replic',B,'burnin',20,'wild');

[T,m]=size(DD);

fprintf('\n Estimating bootstrap sample ... '); kk=1;

IRb = [];

for b=1:B

if b==floor(kk*B/10)

fprintf('%1.0f ',kk);

kk=kk+1;

end

beta = bhat(:,b);

BET = reshape(beta,m,[]);

Sb = Omg(:,:,b);

F = [BET(:,2:end);...

[kron(eye(p-1,p-1),eye(m,m)), zeros(m*(p-1),m)]];

eigen = eig(F);

if ~any(abs(eigen)>1) % is stable?

LR = inv(eye(m)-BET(:,2:m+1)-BET(:,m+2:2*m+1));

% Long-run restriction

RA(m^2,m*(m+1)/2) = -LR(4,4)/LR(4,5);

gammaAb = mlSVAR(RA,gamma0,Sb,T,'tau',.3, ...

'errotol_score',1e-4,'errortol_grad',1e-6,'quiet');

Ab = reshape(RA*gammaAb,m,m);

IRb(:,:,:,b) = VARimpulse(F,Ab,K);

end

end

fprintf(', done!\n');

save SVAR_irf IR IRb DD var_names B

else

load SVAR_irf

end

%% [VII] Plotting IRFs based on structural identification

Kmax = 46;

x = (1:Kmax)';

ll=0.16;

uu=0.84;

sh = [5 4]; % i / Ds

resp = [5 4 2 1]; % i/Ds/pi/y

for ss=1:numel(sh)

for vv=1:numel(resp)

posit = 'NorthEast';

tempb=squeeze(IRb(resp(vv),:,sh(ss),:));

% normalizing

aux = ones(K+1,1)*squeeze(IRb(sh(ss),1,sh(ss),:))';

tempb=tempb./aux;

irf = IR(resp(vv),:,sh(ss))'/IR(sh(ss),1,sh(ss));

if resp(vv)~=4;

irfb = quantile(tempb,.5,2);

irfb_l = quantile(tempb,ll,2);

irfb_u = quantile(tempb,uu,2);

else

irfb = quantile(cumsum(tempb),.5,2);

irfb_l = quantile(cumsum(tempb),ll,2);

irfb_u = quantile(cumsum(tempb),uu,2);

end

irf = irf(x);

irfb=irfb(x); irfb_l=irfb_l(x); irfb_u=irfb_u(x);

miny=min(irfb_l);

maxy=max(irfb_u);

figure(cc); hold on;

dif = irfb_u-irfb_l;

a = area([irfb_l dif],-15);

h=plot(x, irfb);

plot([1 K],[0 0],':k');

set(a(1), 'EdgeColor', 'none', 'FaceColor', 'none');

set(a(2), 'EdgeColor', 'none', 'FaceColor', [0.65 0.65 1]);

set(h, 'Linewidth', 2.5, 'Color', [0 0 0.5]);

% axis

axs=gca;

ylim([min(miny,0)*1.2 maxy*1.2]);

set(axs,'xtick',1:5:Kmax,'Xticklabel',0:5:Kmax-1);

xlabel('Months','fontname','Times','FontWeight', 'normal', 'Fontsize', 22);

ylabel('Percent','fontname','Times','FontWeight', 'normal', 'Fontsize', 22);

formataxis(axs);

xlim([0 Kmax+1]);

if resp(vv)==5

% legend only for one subplot

legend([h;a(2)],'IRF',['Percentiles [' num2str(ll) ' - ' num2str(uu) ']']);

formatlegend(posit);

end

% imprime(['IRF' cod_names{resp(vv)} '_' cod_names{sh(ss)}]);

% imprpdf(['IRF' cod_names{resp(vv)} '_' cod_names{sh(ss)}]);

cc=cc+1;

end

end

toc

%% Extra [Model selection table]

fprintf('\n');

for ii=1:size(ic,1)

fprintf('$%2.0f$ & $%2.4f$ & $%2.4f$ & $%2.4f$ & $%2.4f$ \\\\ \n',ii,ic(ii,:))

end