Use python to develop DSGE II

1451

收藏 2021-04-04

<div class="pinggu_markdown">
<div class="pinggu_markdown__html"><blockquote>
do it best , economy and management. 
Author : Daniel tulips liu . 
University : Renmin Universtiy of China. 
Copyright © Andrea Gazzani (based on the original codes of FORNI, GAMBETTI, LIPPI AND SALA)
</blockquote>
<h2 id="a-dsge-model-use-python">A DSGE Model use Python</h2>
我了显示的格式更美观，必须要打字一部分，占用帖子的前面部分。 
很久没有能发 MARKDOWN 格式的帖子了， 
没想到今天又能发；
挺不错的。
嘿嘿。
NEWS AND NOISE SHOCKS IN HOUSING MARKET
<h2 id="codes">codes</h2>
<pre class=" language-c"><code class="prism language-c">%% Baseline results in "NEWS AND NOISE SHOCKS IN HOUSING MARKET" by Andrea Gazzani (based on the original codes of
%% FORNI, GAMBETTI, LIPPI AND SALA)

clc; close all; clear all;
tic
set(0,'defaultUicontrolFontName','Palatina');set(0,'defaultUitableFontName','Palatina');
set(0,'defaultAxesFontName','Palatina'); set(0,'defaultTextFontName','Palatina');
set(0,'defaultUipanelFontName','Palatina'); set(0,'defaultlinelinewidth',2);

addpath(genpath('Auxiliaries'))

%% Load and prepare data

VARvnames_long = {'Rents';'Housing Prices';'GDP';'Residential Investment';...
'Bba Corp Bond Yield';'S&P Composite Index';}; % define var names
VARvnames = {'rents_fred_raw_rentonly';'hp_av_census';'rgdp'; 'resinv'; 'BAA_yield'; 'sp'};
VARnvar = length(VARvnames);

[xlsdata, xlstext] = xlsread('\DataSet_NewsNoise.xlsx'); % load data
data = Num2NaN(xlsdata(65:end-1,1:end));
dates = xlstext(70:end-1,1);
vnames = xlstext(3,2:end);
for ii=1:length(vnames)
DATA.(vnames{ii}) = data(:,ii);
end

DATA.rents_fred_raw_rentonly = log(100*DATA.rents_fred_raw_rentonly./(DATA.pce_ipd.*DATA.pop)); % real rents
DATA.rents_cpi = log(100*DATA.rents_cpi./(DATA.pce_ipd)); % real rents
DATA.shiller_hp = log(100*DATA.shiller_hp); % real HP - Shiller
DATA.hp_med_census = log(100*DATA.hp_med_census./DATA.ipd); % real HP - Census
DATA.hp_av_census = log(100*DATA.hp_av_census./DATA.ipd); % real HP - Census
DATA.sp = log(100*DATA.sp./DATA.ipd); % real SP
DATA.rgdp = log(100*DATA.rgdp./DATA.pop); % real per capita GDP
DATA.resinv = log(100*DATA.resinv./(DATA.ipd_resinv.*DATA.pop)); % real per capita Res Inv
DATA.hstarts = log(100*DATA.hstarts);
DATA.hpermits = log(100*DATA.hpermits);

nobs = size(data,1);
ENDO = nan(nobs,VARnvar);
for ii=1:VARnvar
ENDO(:,ii) = DATA.(VARvnames{ii});
end

Data = ENDO;

%% Parameters

ll= 40;
npc=0;
nrepli=2000;
maxlags=8;
col2 = [.65 .65 .65];
col1 = [.85 .85 .85];

%storeFEV_surprise = varu./varseries; storeFEV_signal = varalpha./varseries;

varnews0 = squeeze(cumsum(irfs(:,1,:).^2,3)); % variance by news
varnoise0 = squeeze(cumsum(irfs(:,2,:).^2,3)); % variance by noise

varnewsnoise = varnews0 + varnoise0; % news + noise
totvar = totvarnonstruc - varualfa + varnewsnoise; % total var
varnews1 = varnews0./totvar; % variance explained by news
varnoise1 = varnoise0./totvar; % variance explained by noise

tab7= varnews1(varseries,[1 5 9 17 41]);
tab8= varnoise1(varseries,[1 5 9 17 41]);
tab9= varnews1(varseries,[1 5 9 17 41])+varnoise1(varseries,[1 5 9 17 41]);

matrix2lyx(tab7*100, 'tableNEWS.lyx', '%.1f');
matrix2lyx(tab8*100, 'tableNOISE.lyx', '%.1f');
matrix2lyx(tab9*100, 'tableNEWS+NOISE.lyx', '%.1f');
matrix2lyx(tab5*100, 'tableLEARNING.lyx', '%.1f');
matrix2lyx(tab6*100, 'tableSIGNAL.lyx', '%.1f');

%% Bootstrap standard errors

X = Data(:,varseries); % load data

[redAR,sigma,sidui,R2,constant,Yhat,ex]=varestimy(lags,X,1); % estimate VAR
ciao = chol(cov(sidui));
[T N] = size(X);
AA = [];
for t=2:lags+1
AA = [AA squeeze(-redAR(:,:,t))];
end

hor = 40;
jota = [AA ; [eye(N*(lags-1),N*(lags-1)) zeros(N*(lags-1),N)]];
for s = 1:hor
IRF_j = jota^(s-1);
IRF_lev(:,:,s) = IRF_j(1:N,1:N,:)*ciao';
end
temp = sidui*inv(ciao);

standard_shocks = temp;

p = lags;
[A,SIGMA,U,V]=olsvarc(X,p); % estimate VAR in levels

h = 40;
[CI,CI_chol]=boot_luca_MINE2(A,U,X,V,p,h,nrepli,varseries,ind1,ind); % bootstrap

irfold = prctile(CI_chol,[50],4);
irf_ls = irfold;

%% Cumulative FEVD

totvarnonstruc = squeeze(sum(cumsum(irfold.^2,3),2)); % total variance

varu = squeeze(cumsum(irf_ls(:,ind1,:).^2,3)); % variance by learning
varalfa = squeeze(cumsum(irf_ls(:,ind,:).^2,3)); % variance by signal
varualfa = varu + varalfa;

varsignal = varalfa./totvarnonstruc; varsurprise = varu./totvarnonstruc;

irf_nn = prctile(CI,[50],4);

varnews0 = squeeze(cumsum(irf_nn(:,1,:).^2,3)); % variance by news
varnoise0 = squeeze(cumsum(irf_nn(:,2,:).^2,3)); % variance by noise

varnewsnoise = varnews0 + varnoise0; % news + noise
totvar = totvarnonstruc - varualfa + varnewsnoise; % total var
varnews1 = varnews0./totvar; % variance explained by news
varnoise1 = varnoise0./totvar; % variance explained by noise

varnews1(1,1) = 0; varnoise1(1,1) = 0; % implied by theoretical definitions in FGLS
tab_varSIGNAL_CumSum = varsignal(:,[1 5 9 17 41]); tab_varSURPRISE_CumSum = varsurprise(:,[1 5 9 17 41]);
tab_varNEWS_CumSum = varnews1(:,[1 5 9 17 41]); tab_varNOISE_CumSum=varnoise1(:,[1 5 9 17 41]);

ff=1;
figure(ff);
FigSize(30,21)
for t=1:size(varnews1,1)
 subplot(size(varnews1,1)/2,2,t);
 a=bar([(100*varnews1(t,1:ll))' 100*varnoise1(t,1:ll)'],'stacked');grid('on'); %axis tight;
 set(a(1),'facecolor','b');
 set(a(2),'facecolor','r');
 tit(t)=title(ti{t});
 xlab(t)=xlabel('Horizon','FontSize',12);
 ylab(t)=ylabel('% FEV');
 axis tight;
 %axis([0 ll-1 0 100])
 optfont = FigFontOption_Typewriter; optfont.title_weight = 'bold'; FigFont(optfont)
end
leg=legend('News','Noise');
set(ylab,'FontSize',12); set(leg,'FontSize',14); set(tit,'FontSize',14); set(xlab,'FontSize',12);
spaceplots([0.01 0.01 0.01 0.01], [.02 .02]);
set(gcf, 'Color', 'w');

varnews0_p = squeeze(irf_nn(:,1,:).^2); % variance by news
varnoise0_p = squeeze(irf_nn(:,2,:).^2); % variance by noise

varnewsnoise_p = varnews0_p + varnoise0_p; % news + noise
totvar_p = totvarnonstruc_p - varualfa_p + varnewsnoise_p; % total var
varnews1_p = varnews0_p./totvar_p; % variance explained by news
varnoise1_p = varnoise0_p./totvar_p; % variance explained by noise

varlearning_p= varu_p./totvarnonstruc_p; % variance explained by learnig instantaneously
varsignal_p= varalfa_p./totvarnonstruc_p; % variance explained by signal instantaneously, 1 year, 2 year,...
varnews1_p(1,1) = 0; varnoise1_p(1,1) = 0; % implied by theoretical definitions of FGLS

ff=1;
figure(ff);
FigSize(30,21)
for t=1:size(varnews1_p,1)
 subplot(size(varnews1_p,1)/2,2,t);
 a=bar([(100*varnews1_p(t,1:ll))' 100*varnoise1_p(t,1:ll)'],'stacked');grid('on'); %axis tight;
 set(a(1),'facecolor','b');
 set(a(2),'facecolor','r');
 tit(t)=title(ti{t});
 xlab(t)=xlabel('Horizon','FontSize',12);
 ylab(t)=ylabel('% FEV');
 axis tight;
 %axis([0 ll-1 0 100])
 optfont = FigFontOption_Typewriter; optfont.title_weight = 'bold'; FigFont(optfont)
end
leg=legend('News','Noise');
set(ylab,'FontSize',12); set(leg,'FontSize',14); set(tit,'FontSize',14); set(xlab,'FontSize',12);
spaceplots([0.01 0.01 0.01 0.01], [.02 .02]);
set(gcf, 'Color', 'w');

% Store Tables

tab_varSIGNAL = varsignal_p(:,[1 5 9 17 41]); tab_varSURPRISE = varlearning_p(:,[1 5 9 17 41]);
tab_varNEWS = varnews1_p(:,[1 5 9 17 41]); tab_varNOISE=varnoise1_p(:,[1 5 9 17 41]);

matrix2lyx(tab_varNEWS*100, 'NEWS_HbyH.lyx', '%.1f');
matrix2lyx(tab_varNOISE*100, 'NOISE_HbyH.lyx', '%.1f');
matrix2lyx((tab_varNOISE+tab_varNEWS)*100, 'NEWS+NOISE_HbyH.lyx', '%.1f');
matrix2lyx(tab_varSURPRISE*100, 'SURPRISE_HbyH.lyx', '%.1f');
matrix2lyx(tab_varSIGNAL*100, 'SIGNAL_HbyH.lyx', '%.1f');

%% Second ordering

varseries2 = [3 4 5 1 2 6];

% Define indexes for variables

ind1 = 4; ind = 5;

codeData=ones(1,6);codeData(varseries2(:)) = 0;
nr = 10;

% Computation-identification

[irfs2,irfsb2,irf2,irfb2,sigma_a2,sigma_e2,sigma_aboot2,bidielle2,sh2,ssh2,wr2,shocks2] = FAVARNewsNoiseChol_allshocks(Data,...
[Data(:,varseries2) pc(:,1:npc)],codeData,lags,ll,nr,[varseries2 15:15+npc-1],ind1,ind);

[r2,p2,rlo2,rup2]=corrcoef([ssh2 shocks2]);

% Variance Decomposition

irfold2 = irf2;
totvarnonstruc2 = squeeze(sum(cumsum(irfold2.^2,3),2)); % total variance

varu2= squeeze(cumsum(irf2(:,ind1,:).^2,3)); % variance by learning
varalfa2 = squeeze(cumsum(irf2(:,ind,:).^2,3)); % variance by signal
varualfa2 = varu2 + varalfa2;

tab52= varu2(varseries2,[1 5 9 17 41])./totvarnonstruc2(varseries2,[1 5 9 25 41]); % variance explained by learnig instantaneously, 1 year, 2 year,...
tab62= varalfa2(varseries2,[1 5 9 17 41])./totvarnonstruc2(varseries2,[1 5 9 25 41]); % variance explained by signal instantaneously, 1 year, 2 year,...

varnews02 = squeeze(cumsum(irfs2(:,1,:).^2,3)); % variance by news
varnoise02 = squeeze(cumsum(irfs2(:,2,:).^2,3)); % variance by noise

varnewsnoise2 = varnews02 + varnoise02; % news + noise
totvar2 = totvarnonstruc2 - varualfa2 + varnewsnoise2; % total var
varnews12= varnews02./totvar2; % variance explained by news
varnoise12 = varnoise02./totvar2; % variance explained by noise

tab72= varnews12(varseries2,[1 5 9 17 41 ll+1]);
tab82= varnoise12(varseries2,[1 5 9 17 41 ll+1]);
tab92 = varnews12(varseries2,[1 5 9 17 41 ll+1])+varnoise12(varseries2,[1 5 9 17 41 ll+1]);

matrix2lyx(tab72*100, 'NEWS_2ord.lyx', '%.1f');
matrix2lyx(tab82*100, 'NOISE_2ord.lyx', '%.1f');
matrix2lyx((tab92)*100, 'NEWS+NOISE_2ord.lyx', '%.1f');
matrix2lyx(tab52*100, 'SURPRISE_2ord.lyx', '%.1f');
matrix2lyx(tab62*100, 'SIGNAL_2ord.lyx', '%.1f');

ind1 = 1;
ind = 2;
indGDP =3;
indC= 7;
indI= 4;
indFFR=5;
indpi=6;
indD = 8;

ff=2;
DoFigures_Hist

</code></pre>
<h1 id="enjoy">Enjoy</h1>
Yours Daniel Tulips Liu
Renmin University of China , Beijing.
</div>
</div>