adding occbin original files from 20140630

time-shift
Michel Juillard 2015-05-31 16:48:20 +02:00
parent bbfc650f7f
commit c782ca3686
21 changed files with 1745 additions and 0 deletions

View File

@ -0,0 +1,21 @@
% Solve model, generate model IRFs
[zdatalinear zdatapiecewise zdatass oobase_ Mbase_ ] = ...
solve_one_constraint(modnam,modnamstar,...
constraint, constraint_relax,...
shockssequence,irfshock,nperiods,maxiter);
% unpack the IRFs
for i=1:Mbase_.endo_nbr
eval([deblank(Mbase_.endo_names(i,:)),'_uncdifference=zdatalinear(:,i);']);
eval([deblank(Mbase_.endo_names(i,:)),'_difference=zdatapiecewise(:,i);']);
eval([deblank(Mbase_.endo_names(i,:)),'_ss=zdatass(i);']);
end
nparams = size(Mbase_.param_names,1);
for i = 1:nparams
eval([Mbase_.param_names(i,:),'= Mbase_.params(i);']);
end

View File

@ -0,0 +1,21 @@
[zdatalinear zdatapiecewise zdatass oobase_ Mbase_] = solve_two_constraints(...
modnam_00,modnam_10,modnam_01,modnam_11,...
constraint1, constraint2,...
constraint_relax1, constraint_relax2,...
scalefactormod,irfshock,nperiods,curb_retrench,maxiter);
for i=1:Mbase_.endo_nbr
eval([deblank(Mbase_.endo_names(i,:)),'_uncdifference=zdatalinear(:,i);']);
eval([deblank(Mbase_.endo_names(i,:)),'_difference=zdatapiecewise(:,i);']);
eval([deblank(Mbase_.endo_names(i,:)),'_ss=zdatass(i);']);
end
constraint1_difference = process_constraint(constraint1,'_difference',Mbase_.endo_names,0);
constraint2_difference = process_constraint(constraint2,'_difference',Mbase_.endo_names,0);
nparams = size(Mbase_.param_names,1);
for i = 1:nparams
eval([Mbase_.param_names(i,:),'= Mbase_.params(i);']);
end

83
matlab/occbin/get_deriv.m Executable file
View File

@ -0,0 +1,83 @@
function [hm1,h,hl1,j,resid] = get_deriv(M_,ys_)
iy_ = M_.lead_lag_incidence;
it_ = 1;
x = zeros(1,M_.exo_nbr);
% For most models, there are leads, lags and current values of variables
if size(iy_,1)==3
% find non-zero columns of hm1
lag_cols = find(iy_(1,:)~=0);
% find non-zero columns of h
con_cols = find(iy_(2,:));
% find non-zero columns of hl1
lea_cols = find(iy_(3,:));
% If models either lacks leads or lags, iy_ will have two rows
% In this case, we guess that the row with more nonzeros is the row with current variables
elseif size(iy_,1)==2
% if first row has more nonzero entries than the second, assume model lacks lagged variables
if length(find(iy_(1,:)))>length(find(iy_(2,:)))
warning('Model does not have lagged endogenous variables')
con_cols = find(iy_(1,:));
lea_cols = find(iy_(2,:));
lag_cols = [];
else
warning('Model does not have expected future endogenous variables')
lag_cols = find(iy_(1,:));
con_cols = find(iy_(2,:));
lea_cols = [];
end
end
% find number of entries for y vector
ny = length(find(iy_~=0));
% build steady state y
y = ys_(lag_cols);
y = [y;ys_(con_cols)];
y = [y;ys_(lea_cols)];
if ismac
eval(['[resid,g1]=',M_.fname,'_dynamic(y,x, M_.params, ys_, it_);']);
% Older versions of DYNARE for Mac did not include ys_ in the call structure
%eval(['[resid,g1]=',M_.fname,'_dynamic(y,x, M_.params, it_);']);
else
eval(['[resid,g1]=',M_.fname,'_dynamic(y,x, M_.params, ys_, it_);']);
end
hm1=zeros(M_.endo_nbr);
h = hm1;
hl1 = hm1;
j = zeros(M_.endo_nbr,M_.exo_nbr);
% build hm1
nlag_cols = length(lag_cols);
for i=1:nlag_cols
hm1(:,lag_cols(i)) = g1(:,i);
end
% build h
ncon_cols = length(con_cols);
for i=1:ncon_cols
h(:,con_cols(i)) = g1(:,i+nlag_cols);
end
% build hl1
nlea_cols = length(lea_cols);
for i=1:nlea_cols
hl1(:,lea_cols(i)) = g1(:,i+nlag_cols+ncon_cols);
end
for i = 1:M_.exo_nbr;
j(:,i) =g1(:,i+ny);
end

28
matlab/occbin/get_pq.m Executable file
View File

@ -0,0 +1,28 @@
function [p,q]=get_pq(dr_,nstatic,nfwrd);
nvars = size(dr_.ghx,1);
nshocks = size(dr_.ghu,2);
statevar_pos = (nstatic +1):(nvars-nfwrd);
p = zeros(nvars);
% interlace matrix
nnotzero = length(statevar_pos);
for i=1:nnotzero
p(:,statevar_pos(i)) = dr_.ghx(:,i);
end
% reorder p matrix according to order in lgy_
inverse_order = zeros(nvars,1);
for i=1:nvars
inverse_order(i) = find(i==dr_.order_var);
end
p_reordered = zeros(nvars);
q = zeros(nvars,nshocks);
for i=1:nvars
for j=1:nvars
p_reordered(i,j)=p(inverse_order(i),inverse_order(j));
end
q(i,:)=dr_.ghu(inverse_order(i),:);
end
p=p_reordered;

81
matlab/occbin/makechart.m Executable file
View File

@ -0,0 +1,81 @@
function makechart(titlelist,legendlist,figlabel,ylabels,zdata1,zdata2,zdata3)
figure
titlelist = char(strrep(cellstr(titlelist),'_','.'));
ndsets=3; % default, changed below as applicable
if nargin==5
zdata2=nan*zdata1;
zdata3=nan*zdata1;
ndsets =1;
elseif nargin == 6
zdata3 =nan*zdata1;
ndsets=2;
elseif ((nargin>8) | (nargin <=4))
error ('makechart takes 5 to 6 arguments')
end
nobs = size(zdata1,1);
xvalues = (1:nobs)';
nvars = size(titlelist,1);
if nvars==1
nrows=1;
ncols = 1;
elseif nvars==2
nrows =2;
ncols = 1;
elseif (nvars == 3 | nvars ==4)
nrows = 2;
ncols =2;
elseif (nvars==5 |nvars ==6)
nrows = 3;
ncols = 2;
elseif (nvars==7 | nvars==8)
nrows = 4;
ncols = 2;
elseif (nvars==9 | nvars==10)
nrows = 5;
ncols = 2;
else
error('too many variables (makechart)')
end
for i = 1:nvars
subplot(nrows,ncols,i)
h1=plot(xvalues,zdata1(:,i),'b-','linewidth',2); hold on
h1=plot(xvalues,zdata2(:,i),'r--','linewidth',2); hold on
h2=plot(xvalues,zdata3(:,i),'b-','LineWidth',3);
[x0 x1 y10 y11] = pickaxes(xvalues,zdata1(:,i));
[x0 x1 y20 y21] = pickaxes(xvalues,zdata2(:,i));
[x0 x1 y30 y31] = pickaxes(xvalues,zdata3(:,i));
y0 = min([y10,y20,y30]);
y1 = max([y11,y21,y31]);
if y0==y1
y1=y0+1;
end
axis([x0 x1 y0 y1])
set(h1);
if i==1 && isempty(legendlist)==0
legend(legendlist)
text('String',figlabel,'Units','normalized','Position',[1.2 1.24],...
'FontSize',14,'FontWeight','bold','HorizontalAlignment','center');
end
if i==nvars | i==nvars-1
xlabel('Time');
end
% set(gca,'XTick',xtick)
% set(gca,'XTickLabel',xticklabel)
title([num2str(i),'. ',titlelist(i,:)]);
ylabel(ylabels(i,:))
end
% sets printing preferences
%printpref

136
matlab/occbin/makechart9.m Executable file
View File

@ -0,0 +1,136 @@
function makechart9(titlelist,legendlist,figlabel,yearshock,ylabels,...
zdata1,zdata2,zdata3,zdata4,zdata5,zdata6,zdata7)
figure
titlelist = char(strrep(cellstr(titlelist),'_','.'));
ndsets=7; % default, changed below as applicable
if nargin==6
zdata2=nan*zdata1;
zdata3=nan*zdata1;
zdata4=nan*zdata1;
zdata5=nan*zdata1;
zdata6=nan*zdata1;
zdata7=nan*zdata1;
ndsets =1;
elseif nargin==7
zdata3=nan*zdata1;
zdata4=nan*zdata1;
zdata5=nan*zdata1;
zdata6=nan*zdata1;
zdata7=nan*zdata1;
ndsets =2;
elseif nargin == 8
zdata4 =nan*zdata1;
zdata5 =nan*zdata1;
zdata6=nan*zdata1;
zdata7=nan*zdata1;
ndsets=3;
elseif nargin == 9
zdata5 =nan*zdata1;
zdata6=nan*zdata1;
zdata7=nan*zdata1;
ndsets=4;
elseif nargin == 10
zdata6 =nan*zdata1;
zdata7=nan*zdata1;
ndsets=5;
elseif nargin == 11
zdata7=nan*zdata1;
ndsets=6;
elseif ((nargin>=13) | (nargin <=3))
error ('makechart takes 4 to 10 arguments')
end
nobs = size(zdata1,1);
if yearshock>-100
xvalues = yearshock+(0:nobs-1)'/4; % Matteo plot year on x axis
else
xvalues = (1:nobs)'; % Matteo plot year on x axis
end
nvars = size(titlelist,1);
if nvars==1
nrows=1;
ncols = 1;
elseif nvars==2
nrows =2;
ncols = 1;
elseif nvars == 3
nrows = 3;
ncols = 1;
elseif nvars==4
nrows = 2;
ncols = 2;
elseif (nvars==5 | nvars ==6)
nrows = 3;
ncols = 2;
elseif (nvars==7 | nvars==8)
nrows = 4;
ncols = 2;
elseif nvars>8 & nvars<=12;
nrows = 3;
ncols = 4;
elseif nvars>12 & nvars<=15;
nrows = 5;
ncols = 3;
else
error('too many variables (makechart)')
end
for i = 1:nvars
subplot(nrows,ncols,i)
h1=plot(xvalues,zdata1(:,i),'k',...
xvalues,zdata2(:,i),'r',...
xvalues,zdata3(:,i),'b',...
xvalues,zdata4(:,i),'g',...
xvalues,zdata5(:,i),'g',...
xvalues,zdata6(:,i),'c',...
xvalues,zdata7(:,i),'y');
[x0 x1 y10 y11] = pickaxes(xvalues,zdata1(:,i));
[x0 x1 y20 y21] = pickaxes(xvalues,zdata2(:,i));
[x0 x1 y30 y31] = pickaxes(xvalues,zdata3(:,i));
[x0 x1 y40 y41] = pickaxes(xvalues,zdata4(:,i));
[x0 x1 y50 y51] = pickaxes(xvalues,zdata5(:,i));
[x0 x1 y60 y61] = pickaxes(xvalues,zdata6(:,i));
[x0 x1 y70 y71] = pickaxes(xvalues,zdata7(:,i));
grid on
y0 = min([y10,y20,y30,y40,y50,y60,y70]);
y1 = max([y11,y21,y31,y41,y51,y61,y71]);
if y0==y1
y1=y0+1;
end
axis([x0 x1 y0 y1])
set(h1,'linewidth',2);
if i==1
if numel(strvcat(legendlist(1,:)))
h=legend(legendlist,'Location','Northwest');
set(h,'Fontsize',8)
end
end
if i==1
if nvars>3
text('String',figlabel,'Units','normalized','Position',[1.2 1.21],...
'FontSize',13,'FontWeight','bold','HorizontalAlignment','center');
else
text('String',figlabel,'Units','normalized','Position',[0.4 1.24],...
'FontSize',13,'FontWeight','bold','HorizontalAlignment','center');
end
end
%set(gca,'XTick',xtick)
%set(gca,'XTickLabel',xticklabel)
title(titlelist(i,:),'Fontsize',11);
ylabel(ylabels(i,:))
end
% sets printing preferences
%printpref

25
matlab/occbin/map_regime.m Executable file
View File

@ -0,0 +1,25 @@
function [regime regimestart]=map_regimes(violvecbool)
nperiods = length(violvecbool)-1;
% analyse violvec and isolate contiguous periods in the other regime.
regime(1) = violvecbool(1);
regimeindx = 1;
regimestart(1) = 1;
for i=2:nperiods
if violvecbool(i)~=regime(regimeindx)
regimeindx=regimeindx+1;
regime(regimeindx) = violvecbool(i);
regimestart(regimeindx)=i;
end
end
if (regime(1) == 1 & length(regimestart)==1)
warning('Increase nperiods');
end
if (regime(end)==1)
warning('Increase nperiods');
end

63
matlab/occbin/mkdata.m Executable file
View File

@ -0,0 +1,63 @@
function [zdata]=mkdata(nperiods,decrulea,decruleb,endog_,exog_,wishlist,irfshock,scalefactormod,init)
%[nsim, ksim, ysim, isim, csim] = mkdata(nperiods,cofb,endog_)
% given decision rule
neqs = size(endog_,1);
if nargin<9
init = zeros(neqs,1);
end
if nargin<8
scalefactormod=1;
end
if nargin<7
error('Not enough inputs')
end
history = zeros(neqs,nperiods+1);
nshocks = size(irfshock,1);
for i = 1:nshocks
shockpos = strmatch(irfshock(i,:),exog_,'exact');
if ~isempty(shockpos)
irfshockpos(i) = shockpos;
else
error(['Shock ',irfshock(i,:),' is not in the model']);
end
end
% generate data
% history will contain data, the state vector at each period in time will
% be stored columnwise.
history = zeros(neqs,nperiods);
history(:,1)= init;
lengthshock = size(scalefactormod,1);
errvec = zeros(size(exog_,1),1);
for i = 2:nperiods+1
if i<=(lengthshock+1)
for j = 1:nshocks
errvec(irfshockpos(j)) = scalefactormod(i-1,j);
end
history(:,i) = decrulea * history(:,i-1)+decruleb*errvec;
else
% update endogenous variables
history(:,i) = decrulea * history(:,i-1);
end
end
% extract desired variables
nwish=size(wishlist,1);
wishpos = zeros(nwish,1);
history=history';
for i=1:nwish
wishpos(i) = strmatch(wishlist(i,:),endog_,'exact');
end
zdata = history(2:end,wishpos);

View File

@ -0,0 +1,125 @@
function [zdata]=mkdatap_anticipated(nperiods,decrulea,decruleb,...
cof,Jbarmat,cofstar,Jstarbarmat,Dstarbarmat,...
regime,regimestart,violvecbool,...
endog_,exog_,irfshock,scalefactormod,init)
nvars = size(endog_,1);
if nargin<16
init=zeros(nvars,1);
end
if nargin<15;
scalefactormod=1;
end
nshocks = size(irfshock,1);
for i = 1:nshocks
shockpos = strmatch(irfshock(i,:),exog_,'exact');
if ~isempty(shockpos)
irfshockpos(i) = shockpos;
else
error(['Shock ',irfshock(i,:),' is not in the model']);
end
end
nregimes = length(regime);
Cbarmat = cof(:,1:nvars);
Bbarmat = cof(:,nvars+1:2*nvars);
Abarmat = cof(:,2*nvars+1:3*nvars);
% cofstar contains the system for the model when the constraint binds
Cstarbarmat = cofstar(:,1:nvars);
Bstarbarmat = cofstar(:,nvars+1:2*nvars);
Astarbarmat = cofstar(:,2*nvars+1:3*nvars);
% get the time-dependent decision rules
Tmax = regimestart(nregimes)-1; % Tmax is the position of the last period
% when the constraint binds
if Tmax > 0
P = zeros(nvars,nvars,Tmax);
D = zeros(nvars,Tmax);
invmat = inv((Astarbarmat*decrulea+Bstarbarmat));
P(:,:,Tmax) = -invmat*Cstarbarmat;
D(:,Tmax) = -invmat*Dstarbarmat;
% equivalent to pre-multiplying by the inverse above if the target
% matrix is invertible. Otherwise it yields the minimum state solution
%P(:,:,Tmax) = -(Astarbarmat*decrulea+Bstarbarmat)\Cstarbarmat;
%D(:,Tmax) = -(Astarbarmat*decrulea+Bstarbarmat)\Dstarbarmat;
for i = Tmax-1:-1:1
if violvecbool(i)
invmat = inv(Bstarbarmat+Astarbarmat*P(:,:,i+1));
P(:,:,i)=-invmat*Cstarbarmat;
D(:,i) = -invmat*(Astarbarmat*D(:,i+1)+Dstarbarmat);
else
invmat = inv(Bbarmat+Abarmat*P(:,:,i+1));
P(:,:,i)=-invmat*Cbarmat;
D(:,i) = -invmat*(Abarmat*D(:,i+1));
end
end
if Tmax > 1
if violvecbool(1)
E = -invmat*Jstarbarmat;
else
E = -invmat*Jbarmat;
end
else
invmat = inv(Astarbarmat*decrulea+Bstarbarmat);
E = -invmat*Jstarbarmat;
end
end
% generate data
% history will contain data, the state vector at each period in time will
% be stored columnwise.
history = zeros(nvars,nperiods+1);
history(:,1) = init;
errvec = zeros(size(exog_,1),1);
% deal with predetermined conditions
for i = 1:nshocks
errvec(irfshockpos(i)) = scalefactormod(i);
end
% deal with shocks
irfpos =1;
if irfpos <=Tmax
history(:,irfpos+1) = P(:,:,irfpos)* history(:,irfpos)+...
D(:,irfpos) + E*errvec;
else
history(:,irfpos+1) = decrulea*history(:,irfpos)+decruleb*errvec;
end
% all other periods
for irfpos=2:nperiods+1
if irfpos <=Tmax
history(:,irfpos+1) = P(:,:,irfpos)* history(:,irfpos)+...
D(:,irfpos);
else
history(:,irfpos+1) = decrulea*history(:,irfpos);
end
end
history=history';
zdata = history(2:end,:);

View File

@ -0,0 +1,180 @@
function [zdata]=mkdatap_anticipated_2constraints_alt(nperiods,decrulea,decruleb,...
cof,Jbarmat,...
cof10,Jbarmat10,Dbarmat10,...
cof01,Jbarmat01,Dbarmat01,...
cof11,Jbarmat11,Dbarmat11,...
regime1,regimestart1,...
regime2,regimestart2,...
violvecbool,endog_,exog_,...
irfshock,scalefactormod,init)
nvars = size(endog_,1);
if nargin<16
init=zeros(nvars,1);
end
if nargin<15;
scalefactormod=1;
end
nshocks = size(irfshock,1);
for i = 1:nshocks
shockpos = strmatch(irfshock(i,:),exog_,'exact');
if ~isempty(shockpos)
irfshockpos(i) = shockpos;
else
error(['Shock ',irfshock(i,:),' is not in the model']);
end
end
Cbarmat = cof(:,1:nvars);
Bbarmat = cof(:,nvars+1:2*nvars);
Abarmat = cof(:,2*nvars+1:3*nvars);
% cofstar contains the system for the model when the constraint binds
Cbarmat10 = cof10(:,1:nvars);
Bbarmat10 = cof10(:,nvars+1:2*nvars);
Abarmat10 = cof10(:,2*nvars+1:3*nvars);
Cbarmat01 = cof01(:,1:nvars);
Bbarmat01 = cof01(:,nvars+1:2*nvars);
Abarmat01 = cof01(:,2*nvars+1:3*nvars);
Cbarmat11 = cof11(:,1:nvars);
Bbarmat11 = cof11(:,nvars+1:2*nvars);
Abarmat11 = cof11(:,2*nvars+1:3*nvars);
% get the time-dependent decision rules
nregimes1 = length(regime1);
nregimes2 = length(regime2);
Tmax = max([regimestart1(nregimes1) regimestart2(nregimes2)])-1; % Tmax is the position of the last period
% when the constraint binds
if Tmax > 0
P = zeros(nvars,nvars,Tmax);
D = zeros(nvars,Tmax);
% invmat = inv((Astarbarmat*decrulea+Bstarbarmat));
% P(:,:,Tmax) = -invmat*Cstarbarmat;
% D(:,Tmax) = -invmat*Dstarbarmat;
if (violvecbool(Tmax,1) & ~violvecbool(Tmax,2))
%XXX fix next three lines
invmat = inv((Abarmat10*decrulea+Bbarmat10));
P(:,:,Tmax) = -invmat*Cbarmat10;
D(:,Tmax) = -invmat*Dbarmat10;
elseif (violvecbool(Tmax,1) & violvecbool(Tmax,2))
invmat = inv((Abarmat11*decrulea+Bbarmat11));
P(:,:,Tmax) = -invmat*Cbarmat11;
D(:,Tmax) = -invmat*Dbarmat11;
else
invmat = inv((Abarmat01*decrulea+Bbarmat01));
P(:,:,Tmax) = -invmat*Cbarmat01;
D(:,Tmax) = -invmat*Dbarmat01;
end
for i = Tmax-1:-1:1
if (violvecbool(i,1) & ~violvecbool(i,2))
invmat = inv(Bbarmat10+Abarmat10*P(:,:,i+1));
P(:,:,i)=-invmat*Cbarmat10;
D(:,i) = -invmat*(Abarmat10*D(:,i+1)+Dbarmat10);
elseif (~violvecbool(i,1) & violvecbool(i,2))
invmat = inv(Bbarmat01+Abarmat01*P(:,:,i+1));
P(:,:,i)=-invmat*Cbarmat01;
D(:,i) = -invmat*(Abarmat01*D(:,i+1)+Dbarmat01);
elseif (violvecbool(i,1) & violvecbool(i,2))
invmat = inv(Bbarmat11+Abarmat11*P(:,:,i+1));
P(:,:,i)=-invmat*Cbarmat11;
D(:,i) = -invmat*(Abarmat11*D(:,i+1)+Dbarmat11);
else
invmat = inv(Bbarmat+Abarmat*P(:,:,i+1));
P(:,:,i)=-invmat*Cbarmat;
D(:,i) = -invmat*(Abarmat*D(:,i+1));
end
end
% Double check the appropriate invmat in each case
% right now -- inherited from previous loop
if Tmax > 1
if ( ~violvecbool(1,1) & violvecbool(1,2) )
E = -invmat*Jbarmat01;
elseif ( violvecbool(1,1) & ~violvecbool(1,2) )
E = -invmat*Jbarmat10;
elseif ( violvecbool(1,1) & violvecbool(1,2) )
E = -invmat*Jbarmat11;
else
E = -invmat*Jbarmat;
end
else % Tmax is equal to 1
% invmat = inv((Astarbarmat*decrulea+Bstarbarmat));
% E = -invmat*Jstarbarmat;
if ( ~violvecbool(1,1) & violvecbool(1,2) )
invmat = inv((Abarmat01*decrulea+Bbarmat01));
E = -invmat*Jbarmat01;
elseif ( violvecbool(1,1) & violvecbool(1,2) )
invmat = inv((Abarmat11*decrulea+Bbarmat11));
E = -invmat*Jbarmat11;
else
invmat = inv((Abarmat10*decrulea+Bbarmat10));
E = -invmat*Jbarmat10;
end
end
end
% generate data
% history will contain data, the state vector at each period in time will
% be stored columnwise.
history = zeros(nvars,nperiods+1);
history(:,1) = init;
errvec = zeros(size(exog_,1),1);
for i = 1:nshocks
errvec(irfshockpos(i)) = scalefactormod(i);
end
% deal with shocks
irfpos =1;
if irfpos <=Tmax
history(:,irfpos+1) = P(:,:,irfpos)* history(:,irfpos)+...
D(:,irfpos) + E*errvec;
else
history(:,irfpos+1) = decrulea*history(:,irfpos)+decruleb*errvec;
end
% all other periods
for irfpos=2:nperiods+1
if irfpos <=Tmax
history(:,irfpos+1) = P(:,:,irfpos)* history(:,irfpos)+...
D(:,irfpos);
else
history(:,irfpos+1) = decrulea*history(:,irfpos);
end
end
history=history';
zdata = history(2:end,:);

16
matlab/occbin/pickaxes.m Executable file
View File

@ -0,0 +1,16 @@
function [x0,x1,y0,y1] = pickaxes(xvalues,yvalues)
x0=xvalues(1);
nobs = length(xvalues);
x1=xvalues(nobs);
maxy = max(yvalues);
miny = min(yvalues);
y0 = miny - .05*abs(miny);
if (miny>0 & y0<0)
y0 = 0;
end
y1 = maxy + .05*abs(maxy);

View File

@ -0,0 +1,38 @@
% this function looks for occurrences of the endogenous variables in
% endo_names in the input string constraint
% all occurrences of the endogenous variables are appended a suffix
% if the invert_switch is true, the direction of the inequality in the
% constraint is inverted
function constraint1 = process_constraint(constraint,suffix,endo_names,invert_switch)
% create a list of delimiters that can separate parameters and endogenoous
% variables in the string that expresses the constraint
delimiters = char(',',';','(',')','+','-','^','*','/',' ','>','<','=');
% split the string that holds the constraint into tokens
tokens = tokenize(constraint,delimiters);
ntokens = length(tokens);
% search for tokens that match the list of endogenous variables
for i=1:ntokens
if ~isempty(find(strcmp(tokens(i),endo_names)))
% when there is a match with an endogenous variable append the
% suffix
tokens(i) = cellstr([char(tokens(i)),suffix]);
end
% if the invert_switch is true
% reverse the direction of the inequality
if invert_switch
if strcmp(tokens(i),cellstr('>'))
tokens(i) = cellstr('<');
elseif strcmp(tokens(i),cellstr('<'))
tokens(i) = cellstr('>');
end
end
end
% reassemble the tokens to create a string that expresses the constraint
constraint1 = strmerge(tokens);

14
matlab/occbin/setss.m Executable file
View File

@ -0,0 +1,14 @@
% Script that retrieves parameter values once model is solved
nendog = size(Mbase_.endo_names,1);
for i=1:nendog
eval([deblank(Mbase_.endo_names(i,:)) '_ss = oo_.dr.ys(i); ']);
end
nparams = size(Mbase_.param_names);
for i = 1:nparams
eval([Mbase_.param_names(i,:),'= M_.params(i);']);
end

View File

@ -0,0 +1,50 @@
function [zdata_, zdatass_, oobase_, Mbase_ ] = ...
solve_no_constraint(modnam,...
shockssequence,irfshock,nperiods)
global M_ oo_
errlist = [];
% solve model
eval(['dynare ',modnam,' noclearall'])
oobase_ = oo_;
Mbase_ = M_;
ys_ = oobase_.dr.ys;
nvars = numel(ys_);
zdatass_ = ys_ ;
for i=1:Mbase_.endo_nbr
eval([deblank(Mbase_.endo_names(i,:)) '_ss = oo_.dr.ys(i); ']);
end
for i = 1:size(Mbase_.param_names)
eval([Mbase_.param_names(i,:),'= M_.params(i);']);
end
[hm1,h,hl1,Jbarmat] = get_deriv(Mbase_,ys_);
cof = [hm1,h,hl1];
[decrulea,decruleb]=get_pq(oobase_.dr);
endog_ = M_.endo_names;
exog_ = M_.exo_names;
nvars = numel(Mbase_.endo_nbr);
nshocks = size(shockssequence,1);
init = zeros(nvars,1);
wishlist = endog_;
nwishes = size(wishlist,1);
zdata_ = mkdata(nperiods,decrulea,decruleb,endog_,exog_,wishlist,irfshock,shockssequence);

View File

@ -0,0 +1,48 @@
function [zdata oobase_ Mbase_ ] = ...
solve_no_constraint_noclear(modnam,...
shockssequence,irfshock,nperiods);
global M_ oo_
errlist = [];
% solve model
eval(['dynare ',modnam,' nolog']);
oobase_ = oo_;
Mbase_ = M_;
ys_ = oobase_.dr.ys;
for i=1:Mbase_.endo_nbr
eval([deblank(Mbase_.endo_names(i,:)) '_ss = oo_.dr.ys(i); ']);
end
for i = 1:size(Mbase_.param_names)
eval([Mbase_.param_names(i,:),'= M_.params(i);']);
end
setss
[hm1,h,hl1,Jbarmat] = get_deriv(Mbase_,ys_);
cof = [hm1,h,hl1];
[decrulea,decruleb]=get_pq(oobase_.dr);
endog_ = M_.endo_names;
exog_ = M_.exo_names;
nvars = numel(Mbase_.endo_nbr);
nshocks = size(shockssequence,1);
init = zeros(nvars,1);
wishlist = endog_;
nwishes = size(wishlist,1);
zdata = mkdata(nperiods,decrulea,decruleb,endog_,exog_,wishlist,irfshock,shockssequence);

View File

@ -0,0 +1,192 @@
% solve_one_constraint [zdatalinear zdatapiecewise zdatass oo base M base] = solve one constraint(modnam, modnamstar, constraint, constraint relax, shockssequence, irfshock, nperiods, maxiter, init);
%
% Inputs:
% modnam: name of .mod file for the reference regime (excludes the .mod extension).
% modnamstar: name of .mod file for the alternative regime (excludes the .mod exten- sion).
% constraint: the constraint (see notes 1 and 2 below). When the condition in constraint evaluates to true, the solution switches from the reference to the alternative regime.
% constraint relax: when the condition in constraint relax evaluates to true, the solution returns to the reference regime.
% shockssequence: a sequence of unforeseen shocks under which one wants to solve the model (size T×nshocks).
% irfshock: label for innovation for IRFs, from Dynare .mod file (one or more of the ?varexo?).
% nperiods: simulation horizon (can be longer than the sequence of shocks defined in shockssequence; must be long enough to ensure convergence back to the reference model at the end of the simulation horizon and may need to be varied depending on the sequence of shocks).
% maxiter: maximum number of iterations allowed for the solution algorithm (20 if not specified).
% init: the initial position for the vector of state variables, in deviation from steady state (if not specified, the default is steady state). The ordering follows the definition order in the .mod files.
%
% Outputs:
% zdatalinear: an array containing paths for all endogenous variables ignoring the occasionally binding constraint (the linear solution), in deviation from steady state. Each column is a variable, the order is the definition order in the .mod files.
% zdatapiecewise: an array containing paths for all endogenous variables satisfying the occasionally binding constraint (the occbin/piecewise solution), in deviation from steady state. Each column is a variable, the order is the definition order in the .mod files.
% zdatass: theinitialpositionforthevectorofstatevariables,indeviationfromsteady state (if not specified, the default is a vectors of zero implying that the initial conditions coincide with the steady state). The ordering follows the definition order in the .mod files.
% oobase,Mbase: structures produced by Dynare for the reference model ? see Dynare User Guide.
% Log of changes:
% 6/17/2013 -- Luca added a trailing underscore to local variables in an
% attempt to avoid conflicts with parameter names defined in the .mod files
% to be processed.
% 6/17/2013 -- Luca replaced external .m file setss.m
function [zdatalinear_ zdatapiecewise_ zdatass_ oobase_ Mbase_ ] = ...
solve_one_constraint(modnam_,modnamstar_,...
constraint_, constraint_relax_,...
shockssequence_,irfshock_,nperiods_,maxiter_,init_)
global M_ oo_
errlist_ = [];
% solve the reference model linearly
eval(['dynare ',modnam_,' noclearall nolog '])
oobase_ = oo_;
Mbase_ = M_;
% import locally the values of parameters assigned in the reference .mod
% file
for i_indx_ = 1:Mbase_.param_nbr
eval([Mbase_.param_names(i_indx_,:),'= M_.params(i_indx_);']);
end
% Create steady state values of the variables if needed for processing the constraint
for i=1:Mbase_.endo_nbr
eval([deblank(Mbase_.endo_names(i,:)) '_ss = oobase_.dr.ys(i); ']);
end
% parse the .mod file for the alternative regime
eval(['dynare ',modnamstar_,' noclearall nolog '])
oostar_ = oo_;
Mstar_ = M_;
% check inputs
if ~strcmp(Mbase_.endo_names,Mstar_.endo_names)
error('The two .mod files need to have exactly the same endogenous variables declared in the same order')
end
if ~strcmp(Mbase_.exo_names,Mstar_.exo_names)
error('The two .mod files need to have exactly the same exogenous variables declared in the same order')
end
if ~strcmp(Mbase_.param_names,Mstar_.param_names)
warning('The parameter list does not match across .mod files')
end
% ensure that the two models have the same parameters
% use the parameters for the base model.
Mstar_.params = Mbase_.params;
nvars_ = Mbase_.endo_nbr;
zdatass_ = oobase_.dr.ys;
% get the matrices holding the first derivatives for the model
% each regime is treated separately
[hm1_,h_,hl1_,Jbarmat_] = get_deriv(Mbase_,zdatass_);
cof_ = [hm1_,h_,hl1_];
[hm1_,h_,hl1_,Jstarbarmat_,resid_] = get_deriv(Mstar_,zdatass_);
cofstar_ = [hm1_,h_,hl1_];
Dstartbarmat_ = resid_;
[decrulea_,decruleb_]=get_pq(oobase_.dr);
endog_ = M_.endo_names;
exog_ = M_.exo_names;
% processes the constraints specified in the call to this function
% uppend a suffix to each endogenous variable
constraint_difference_ = process_constraint(constraint_,'_difference',Mbase_.endo_names,0);
constraint_relax_difference_ = process_constraint(constraint_relax_,'_difference',Mbase_.endo_names,0);
nshocks_ = size(shockssequence_,1);
% if necessary, set default values for optional arguments
if ~exist('init_')
init_ = zeros(nvars_,1);
end
if ~exist('maxiter_')
maxiter_ = 20;
end
if ~exist('nperiods_')
nperiods_ = 100;
end
% set some initial conditions and loop through the shocks
% period by period
init_orig_ = init_;
zdatapiecewise_ = zeros(nperiods_,nvars_);
wishlist_ = endog_;
nwishes_ = size(wishlist_,1);
violvecbool_ = zeros(nperiods_+1,1);
for ishock_ = 1:nshocks_
changes_=1;
iter_ = 0;
while (changes_ & iter_<maxiter_)
iter_ = iter_ +1;
% analyze when each regime starts based on current guess
[regime regimestart]=map_regime(violvecbool_);
% get the hypothesized piece wise linear solution
[zdatalinear_]=mkdatap_anticipated(nperiods_,decrulea_,decruleb_,...
cof_,Jbarmat_,cofstar_,Jstarbarmat_,Dstartbarmat_,...
regime,regimestart,violvecbool_,...
endog_,exog_,irfshock_,shockssequence_(ishock_,:),init_);
for i_indx_=1:nwishes_
eval([deblank(wishlist_(i_indx_,:)),'_difference=zdatalinear_(:,i_indx_);']);
end
newviolvecbool_ = eval(constraint_difference_);
relaxconstraint_ = eval(constraint_relax_difference_);
% check if changes to the hypothesis of the duration for each
% regime
if (max(newviolvecbool_-violvecbool_>0)) | sum(relaxconstraint_(find(violvecbool_==1))>0)
changes_ = 1;
else
changes_ = 0;
end
violvecbool_ = (violvecbool_|newviolvecbool_)-(relaxconstraint_ & violvecbool_);
end
init_ = zdatalinear_(1,:);
zdatapiecewise_(ishock_,:)=init_;
init_= init_';
% reset violvecbool_ for next period's shock -- this resetting is
% consistent with expecting no additional shocks
violvecbool_=[violvecbool_(2:end);0];
end
% if necessary, fill in the rest of the path with the remainder of the
% last IRF computed.
zdatapiecewise_(ishock_+1:end,:)=zdatalinear_(2:nperiods_-ishock_+1,:);
% get the linear responses
zdatalinear_ = mkdata(max(nperiods_,size(shockssequence_,1)),...
decrulea_,decruleb_,endog_,exog_,...
wishlist_,irfshock_,shockssequence_,init_orig_);
if changes_ ==1
display('Did not converge -- increase maxiter_')
end

View File

@ -0,0 +1,200 @@
% solve_one_constraint [zdatalinear zdatapiecewise zdatass oo base M base] = solve one constraint(modnam, modnamstar, constraint, constraint relax, shockssequence, irfshock, nperiods, maxiter, init);
%
% Inputs:
% modnam: name of .mod file for the reference regime (excludes the .mod extension).
% modnamstar: name of .mod file for the alternative regime (excludes the .mod exten- sion).
% constraint: the constraint (see notes 1 and 2 below). When the condition in constraint evaluates to true, the solution switches from the reference to the alternative regime.
% constraint relax: when the condition in constraint relax evaluates to true, the solution returns to the reference regime.
% shockssequence: a sequence of unforeseen shocks under which one wants to solve the model (size T×nshocks).
% irfshock: label for innovation for IRFs, from Dynare .mod file (one or more of the ?varexo?).
% nperiods: simulation horizon (can be longer than the sequence of shocks defined in shockssequence; must be long enough to ensure convergence back to the reference model at the end of the simulation horizon and may need to be varied depending on the sequence of shocks).
% maxiter: maximum number of iterations allowed for the solution algorithm (20 if not specified).
% init: the initial position for the vector of state variables, in deviation from steady state (if not specified, the default is steady state). The ordering follows the definition order in the .mod files.
%
% Outputs:
% zdatalinear: an array containing paths for all endogenous variables ignoring the occasionally binding constraint (the linear solution), in deviation from steady state. Each column is a variable, the order is the definition order in the .mod files.
% zdatapiecewise: an array containing paths for all endogenous variables satisfying the occasionally binding constraint (the occbin/piecewise solution), in deviation from steady state. Each column is a variable, the order is the definition order in the .mod files.
% zdatass: theinitialpositionforthevectorofstatevariables,indeviationfromsteady state (if not specified, the default is a vectors of zero implying that the initial conditions coincide with the steady state). The ordering follows the definition order in the .mod files.
% oobase,Mbase: structures produced by Dynare for the reference model ? see Dynare User Guide.
% Log of changes:
% 6/17/2013 -- Luca added a trailing underscore to local variables in an
% attempt to avoid conflicts with parameter names defined in the .mod files
% to be processed.
% 6/17/2013 -- Luca replaced external .m file setss.m
function [zdatalinear_ zdatapiecewise_ zdatass_ oobase_ Mbase_ ] = ...
solve_one_constraint(modnam_,modnamstar_,...
constraint_, constraint_relax_,...
shockssequence_,irfshock_,nperiods_,maxiter_,init_)
global M_ oo_
errlist_ = [];
% solve the reference model linearly
eval(['dynare ',modnam_,' noclearall nolog '])
oobase_ = oo_;
Mbase_ = M_;
% import locally the values of parameters assigned in the reference .mod
% file
for i_indx_ = 1:Mbase_.param_nbr
eval([Mbase_.param_names(i_indx_,:),'= M_.params(i_indx_);']);
end
% Create steady state values of the variables if needed for processing the constraint
for i=1:Mbase_.endo_nbr
eval([deblank(Mbase_.endo_names(i,:)) '_ss = oobase_.dr.ys(i); ']);
end
% parse the .mod file for the alternative regime
eval(['dynare ',modnamstar_,' noclearall nolog '])
oostar_ = oo_;
Mstar_ = M_;
% check inputs
if ~strcmp(Mbase_.endo_names,Mstar_.endo_names)
error('The two .mod files need to have exactly the same endogenous variables declared in the same order')
end
if ~strcmp(Mbase_.exo_names,Mstar_.exo_names)
error('The two .mod files need to have exactly the same exogenous variables declared in the same order')
end
if ~strcmp(Mbase_.param_names,Mstar_.param_names)
warning('The parameter list does not match across .mod files')
end
% ensure that the two models have the same parameters
% use the parameters for the base model.
Mstar_.params = Mbase_.params;
nvars_ = Mbase_.endo_nbr;
zdatass_ = oobase_.dr.ys;
% get the matrices holding the first derivatives for the model
% each regime is treated separately
[hm1_,h_,hl1_,Jbarmat_] = get_deriv(Mbase_,zdatass_);
cof_ = [hm1_,h_,hl1_];
[hm1_,h_,hl1_,Jstarbarmat_,resid_] = get_deriv(Mstar_,zdatass_);
cofstar_ = [hm1_,h_,hl1_];
Dstartbarmat_ = resid_;
if isfield(Mbase_,'nfwrd')
% the latest Dynare distributions have moved nstatic and nfwrd
[decrulea_,decruleb_]=get_pq(oobase_.dr,Mbase_.nstatic,Mbase_.nfwrd);
else
[decrulea_,decruleb_]=get_pq(oobase_.dr,oobase_.dr.nstatic,oobase_.dr.nfwrd);
end
endog_ = M_.endo_names;
exog_ = M_.exo_names;
% processes the constraints specified in the call to this function
% uppend a suffix to each endogenous variable
constraint_difference_ = process_constraint(constraint_,'_difference',Mbase_.endo_names,0);
constraint_relax_difference_ = process_constraint(constraint_relax_,'_difference',Mbase_.endo_names,0);
nshocks_ = size(shockssequence_,1);
% if necessary, set default values for optional arguments
if ~exist('init_')
init_ = zeros(nvars_,1);
end
if ~exist('maxiter_')
maxiter_ = 20;
end
if ~exist('nperiods_')
nperiods_ = 100;
end
% set some initial conditions and loop through the shocks
% period by period
init_orig_ = init_;
zdatapiecewise_ = zeros(nperiods_,nvars_);
wishlist_ = endog_;
nwishes_ = size(wishlist_,1);
violvecbool_ = zeros(nperiods_+1,1);
for ishock_ = 1:nshocks_
changes_=1;
iter_ = 0;
while (changes_ & iter_<maxiter_)
iter_ = iter_ +1;
% analyze when each regime starts based on current guess
[regime regimestart]=map_regime(violvecbool_);
% get the hypothesized piece wise linear solution
[zdatalinear_]=mkdatap_anticipated(nperiods_,decrulea_,decruleb_,...
cof_,Jbarmat_,cofstar_,Jstarbarmat_,Dstartbarmat_,...
regime,regimestart,violvecbool_,...
endog_,exog_,irfshock_,shockssequence_(ishock_,:),init_);
for i_indx_=1:nwishes_
eval([deblank(wishlist_(i_indx_,:)),'_difference=zdatalinear_(:,i_indx_);']);
end
newviolvecbool_ = eval(constraint_difference_);
relaxconstraint_ = eval(constraint_relax_difference_);
% check if changes to the hypothesis of the duration for each
% regime
if (max(newviolvecbool_-violvecbool_>0)) | sum(relaxconstraint_(find(violvecbool_==1))>0)
changes_ = 1;
else
changes_ = 0;
end
violvecbool_ = (violvecbool_|newviolvecbool_)-(relaxconstraint_ & violvecbool_);
end
init_ = zdatalinear_(1,:);
zdatapiecewise_(ishock_,:)=init_;
init_= init_';
% reset violvecbool_ for next period's shock -- this resetting is
% consistent with expecting no additional shocks
violvecbool_=[violvecbool_(2:end);0];
end
% if necessary, fill in the rest of the path with the remainder of the
% last IRF computed.
zdatapiecewise_(ishock_+1:end,:)=zdatalinear_(2:nperiods_-ishock_+1,:);
% get the linear responses
zdatalinear_ = mkdata(max(nperiods_,size(shockssequence_,1)),...
decrulea_,decruleb_,endog_,exog_,...
wishlist_,irfshock_,shockssequence_,init_orig_);
if changes_ ==1
display('Did not converge -- increase maxiter_')
end

View File

@ -0,0 +1,305 @@
% [zdatalinear zdatapiecewise zdatass oo 00 M 00] = solve two constraints(modnam 00,modnam 10,modnam 01,modnam 11,... constraint1, constraint2,... constraint relax1, constraint relax2,... shockssequence,irfshock,nperiods,curb retrench,maxiter,init);
%
% Inputs:
% modnam 00: name of the .mod file for reference regime (excludes the .mod extension). modnam10: name of the .mod file for the alternative regime governed by the first
% constraint.
% modnam01: name of the .mod file for the alternative regime governed by the second constraint.
% modnam 11: name of the .mod file for the case in which both constraints force a switch to their alternative regimes.
% constraint1: the first constraint (see notes 1 and 2 below). If constraint1 evaluates to true, then the solution switches to the alternative regime for condition 1. In thatcase, if constraint2 (described below) evaluates to false, then the model solution switches to enforcing the conditions for an equilibrium in modnam 10. Otherwise, if constraint2 also evaluates to true, then the model solution switches to enforcing the conditions for an equilibrium in modnam 11.
% constraint relax1: when the condition in constraint relax1 evaluates to true, the solution returns to the reference regime for constraint1.
% constraint2: the second constraint (see notes 1 and 2 below). constraint relax2: when the condition in constraint relax2 evaluates to true, the
% solution returns to the reference regime for constraint2. shockssequence: a sequence of unforeseen shocks under which one wants to solve the
% model
% irfshock: label for innovation for IRFs, from Dynare .mod file (one or more of the ?varexo?)
% nperiods: simulation horizon (can be longer than the sequence of shocks defined in shockssequence; must be long enough to ensure convergence back to the reference model at the end of the simulation horizon and may need to be varied depending on the sequence of shocks).
% curb retrench: a scalar equal to 0 or 1. Default is 0. When set to 0, it updates the guess based of regimes based on the previous iteration. When set to 1, it updates in a manner similar to a Gauss-Jacobi scheme, slowing the iterations down by updating the guess of regimes only one period at a time.
% maxiter: maximum number of iterations allowed for the solution algorithm (20 if not specified).
% init: the initial position for the vector of state variables, in deviation from steady state (if not specified, the default is a vector of zero implying that the initial conditions coincide with the steady state). The ordering follows the definition order in the .mod files.
%
% Outputs:
% zdatalinear: an array containing paths for all endogenous variables ignoring the occasionally binding constraint (the linear solution), in deviation from steady state. Each column is a variable, the order is the definition order in the .mod files.
% zdatapiecewise: an array containing paths for all endogenous variables satisfying the occasionally binding constraint (the occbin/piecewise solution), in deviation from steady state. Each column is a variable, the order is the definition order in the .mod files.
% zdatass: a vector that holds the steady state values of the endogenous variables ( following the definition order in the .mod file).
% oo00 , M00 : structures produced by Dynare for the reference model ? see Dynare User Guide.
% Log of changes
% 6/17/2013 -- Luca added a trailing underscore to local variables in an
% attempt to avoid conflicts with parameter names defined in the .mod files
% to be processed.
% 6/17/2013 -- Luca replaced external .m file setss.m
function [ zdatalinear_ zdatapiecewise_ zdatass_ oo00_ M00_ ] = ...
solve_two_constraints(modnam_00_,modnam_10_,modnam_01_,modnam_11_,...
constrain1_, constrain2_,...
constraint_relax1_, constraint_relax2_,...
shockssequence_,irfshock_,nperiods_,curb_retrench_,maxiter_,init_)
global M_ oo_
% solve model
eval(['dynare ',modnam_00_,' noclearall nolog'])
oo00_ = oo_;
M00_ = M_;
for i=1:M00_.endo_nbr
eval([deblank(M00_.endo_names(i,:)) '_ss = oo00_.dr.ys(i); ']);
end
for i_indx_ = 1:M00_.param_nbr
eval([M00_.param_names(i_indx_,:),'= M00_.params(i_indx_);']);
end
eval(['dynare ',modnam_10_,' noclearall'])
oo10_ = oo_;
M10_ = M_;
eval(['dynare ',modnam_01_,' noclearall'])
oo01_ = oo_;
M01_ = M_;
eval(['dynare ',modnam_11_,' noclearall'])
oo11_ = oo_;
M11_ = M_;
% do some error checking
% check inputs
if ~strcmp(M00_.endo_names,M10_.endo_names)
error([modnam_00_,' and ',modnam_10_,' need to have exactly the same endogenous variables and they need to be declared in the same order'])
end
if ~strcmp(M00_.exo_names,M10_.exo_names)
error([modnam_00_,' and ',modnam_10_,' need to have exactly the same exogenous variables and they need to be declared in the same order'])
end
if ~strcmp(M00_.param_names,M10_.param_names)
warning(['The parameter list does not match across the files ',modnam_00_,' and ',modnam_10_])
end
if ~strcmp(M00_.endo_names,M01_.endo_names)
error([modnam_00,' and ',modnam_01_,' need to have exactly the same endogenous variables and they need to be declared in the same order'])
end
if ~strcmp(M00_.exo_names,M01_.exo_names)
error([modnam_00_,' and ',modnam_01_,' need to have exactly the same exogenous variables and they need to be declared in the same order'])
end
if ~strcmp(M00_.param_names,M01_.param_names)
warning(['The parameter list does not match across the files ',modnam_00_,' and ',modnam_01_])
end
if ~strcmp(M00_.endo_names,M11_.endo_names)
error([modnam_00_,' and ',modnam_11_,' need to have exactly the same endogenous variables and they need to be declared in the same order'])
end
if ~strcmp(M00_.exo_names,M11_.exo_names)
error([modnam_00_,' and ',modnam_11_,' need to have exactly the same exogenous variables and they need to be declared in the same order'])
end
if ~strcmp(M00_.param_names,M11_.param_names)
warning(['The parameter list does not match across the files ',modnam_00_,' and ',modnam_11_])
end
nvars_ = M00_.endo_nbr;
zdatass_ = oo00_.dr.ys;
[hm1_,h_,hl1_,Jbarmat_] = get_deriv(M00_,zdatass_);
cof_ = [hm1_,h_,hl1_];
M10_.params = M00_.params;
[hm1_,h_,hl1_,Jbarmat10_,resid_] = get_deriv(M10_,zdatass_);
cof10_ = [hm1_,h_,hl1_];
Dbarmat10_ = resid_;
M01_.params = M00_.params;
[hm1_,h_,hl1_,Jbarmat01_,resid_] = get_deriv(M01_,zdatass_);
cof01_ = [hm1_,h_,hl1_];
Dbarmat01_ = resid_;
M11_.params = M00_.params;
[hm1_,h_,hl1_,Jbarmat11_,resid_] = get_deriv(M11_,zdatass_);
cof11_ = [hm1_,h_,hl1_];
Dbarmat11_ = resid_;
if isfield(M00_,'nfwrd') % needed for bakward compatibility with older Dynare releases
[decrulea,decruleb]=get_pq(oo00_.dr,M00_.nstatic,M00_.nfwrd);
else
[decrulea,decruleb]=get_pq(oo00_.dr,oo00_.dr.nstatic,oo00_.dr.nfwrd);
end
endog_ = M00_.endo_names;
exog_ = M00_.exo_names;
% processes the constrain so as to uppend a suffix to each
% endogenous variables
constraint1_difference_ = process_constraint(constrain1_,'_difference',M00_.endo_names,0);
% when the last argument in process_constraint is set to 1, the
% direction of the inequality in the constraint is inverted
constraint_relax1_difference_ = process_constraint(constraint_relax1_,'_difference',M00_.endo_names,0);
% processes the constrain so as to uppend a suffix to each
% endogenous variables
constraint2_difference_ = process_constraint(constrain2_,'_difference',M00_.endo_names,0);
% when the last argument in process_constraint is set to 1, the
% direction of the inequality in the constraint is inverted
constraint_relax2_difference_ = process_constraint(constraint_relax2_,'_difference',M00_.endo_names,0);
nshocks = size(shockssequence_,1);
if ~exist('init_')
init_ = zeros(nvars_,1);
end
if ~exist('maxiter_')
maxiter_ = 20;
end
if ~exist('curb_retrench_')
curb_retrench_ = 0;
end
init_orig_ = init_;
zdatapiecewise_ = zeros(nperiods_,nvars_);
violvecbool_ = zeros(nperiods_+1,2); % This sets the first guess for when
% the constraints are going to hold.
% The variable is a boolean with two
% columns. The first column refers to
% constrain1_; the second to
% constrain2_.
% Each row is a period in time.
% If the boolean is true it indicates
% the relevant constraint is expected
% to evaluate to true.
% The default initial guess is
% consistent with the base model always
% holding -- equivalent to the linear
% solution.
wishlist_ = endog_;
nwishes_ = size(wishlist_,1);
for ishock_ = 1:nshocks
changes_=1;
iter_ = 0;
while (changes_ & iter_<maxiter_)
iter_ = iter_ +1;
% analyse violvec and isolate contiguous periods in the other
% regime.
[regime1 regimestart1]=map_regime(violvecbool_(:,1));
[regime2 regimestart2]=map_regime(violvecbool_(:,2));
[zdatalinear_]=mkdatap_anticipated_2constraints(nperiods_,decrulea,decruleb,...
cof_,Jbarmat_,...
cof10_,Jbarmat10_,Dbarmat10_,...
cof01_,Jbarmat01_,Dbarmat01_,...
cof11_,Jbarmat11_,Dbarmat11_,...
regime1,regimestart1,...
regime2,regimestart2,...
violvecbool_,endog_,exog_,...
irfshock_,shockssequence_(ishock_,:),init_);
for i_indx_=1:nwishes_
eval([deblank(wishlist_(i_indx_,:)),'_difference=zdatalinear_(:,i_indx_);']);
end
newviolvecbool1_ = eval(constraint1_difference_);
relaxconstraint1_ = eval(constraint_relax1_difference_);
newviolvecbool2_ = eval(constraint2_difference_);
relaxconstraint2_ = eval(constraint_relax2_difference_);
newviolvecbool_ = [newviolvecbool1_;newviolvecbool2_];
relaxconstraint_ = [relaxconstraint1_;relaxconstraint2_];
% check if changes_
if (max(newviolvecbool_(:)-violvecbool_(:)>0)) | sum(relaxconstraint_(find(violvecbool_==1))>0)
changes_ = 1;
else
changes_ = 0;
end
if curb_retrench_ % apply Gauss-Sidel idea of slowing down the change in the guess
% for the constraint -- only relax one
% period at a time starting from the last
% one when each of the constraints is true.
retrench = 0*violvecbool_(:);
if ~isempty(find(relaxconstraint1_ & violvecbool_(:,1)))
retrenchpos = max(find(relaxconstraint1_ & violvecbool_(:,1)));
retrench(retrenchpos) = 1;
end
if ~isempty(find(relaxconstraint2_ & violvecbool_(:,2)))
retrenchpos = max(find(relaxconstraint2_ & violvecbool_(:,2)));
retrench(retrenchpos+nperiods_+1) = 1;
end
violvecbool_ = (violvecbool_(:) | newviolvecbool_(:))-retrench(:);
else
violvecbool_ = (violvecbool_(:) | newviolvecbool_(:))-(relaxconstraint_(:) & violvecbool_(:));
end
violvecbool_ = reshape(violvecbool_,nperiods_+1,2);
end
if changes_ ==1
display('Did not converge -- increase maxiter')
end
init_ = zdatalinear_(1,:);
zdatapiecewise_(ishock_,:)=init_;
init_= init_';
% update the guess for constraint violations for next period
% update is consistent with expecting no additional shocks next period
violvecbool_=[violvecbool_(2:end,:);zeros(1,2)];
end
zdatapiecewise_(ishock_+1:end,:)=zdatalinear_(2:nperiods_-ishock_+1,:);
zdatalinear_ = mkdata(nperiods_,decrulea,decruleb,endog_,exog_,wishlist_,irfshock_,shockssequence_,init_orig_);

9
matlab/occbin/strmerge.m Executable file
View File

@ -0,0 +1,9 @@
function string = strmerge(tokens)
ntokens = length(tokens);
string = char(tokens(1));
for i=2:ntokens
string = [string,char(tokens(i))];
end

55
matlab/occbin/tokenize.m Executable file
View File

@ -0,0 +1,55 @@
function tokens = tokenize(source,delimiter)
% syntax
% tokens = tokenize(source,delimiters)
%
% source is a string to be broken into tokens
% delimiters is a character array of single character delimiters
% tokens is a cell string array containing the tokens
posdelims = [];
% assumes that delimiter cannot be in the first position or the last
% position
ndelimiters = size(delimiter,1);
for i=1:ndelimiters
newpositions = strfind(source,delimiter(i,:));
if ~isempty(newpositions)
posdelims =[posdelims, newpositions];
end
end
% reorder posdelims in ascending order
posdelims = sort(posdelims);
if isempty(posdelims)
tokens = cellstr(source);
else
ndelims = length(posdelims);
% build positions for substrings
delims = zeros(ndelims+1,2);
for i=1:ndelims+1;
if i==1
if posdelims(1) == 1
tokens = cellstr(source(1));
else
delims(i,:) = [1,posdelims(i)-1];
tokens = cellstr(source([delims(i,1):delims(i,2)]));
tokens = [tokens, source(posdelims(i))];
end
elseif i==ndelims+1
if (posdelims(i-1) < length(source))
delims(i,:) = [posdelims(i-1)+1,length(source)];
tokens = [tokens, cellstr(source([delims(i,1):delims(i,2)]))];
end
else
if posdelims(i)>posdelims(i-1)+1
delims(i,:) = [posdelims(i-1)+1,posdelims(i)-1];
tokens = [tokens, cellstr(source([delims(i,1):delims(i,2)]))];
end
tokens = [tokens, source(posdelims(i))];
end
end
end

55
tests/occbin/dynrbc.mod Executable file
View File

@ -0,0 +1,55 @@
// variables
var a, c, i, k, kprev, lambdak;
// innovations to shock processes
varexo erra;
// parameters
parameters ALPHA, DELTAK, BETA, GAMMAC, RHOA, PHII, PHIK, PSI, PSINEG, ISS, KSS ;
model;
# zkss = ((1/BETA-1+DELTAK)/ALPHA)^(1/(ALPHA-1));
# zcss = -DELTAK*zkss + zkss^ALPHA;
# ziss = DELTAK*zkss;
# zuss = (zcss^(1-GAMMAC)-1)/(1-GAMMAC);
# zvss = zuss/(1-BETA);
/////////////////////////////////////////////////////////////////
// 1.
-exp(c)^(-GAMMAC)*(1+2*PSI*(exp(k)/exp(k(-1))-1)/exp(k(-1)))
+ BETA*exp(c(1))^(-GAMMAC)*((1-DELTAK)-2*PSI*(exp(k(1))/exp(k)-1)*
(-exp(k(1))/exp(k)^2)+ALPHA*exp(a(1))*exp(k)^(ALPHA-1))=
-lambdak+BETA*(1-DELTAK+PHIK)*lambdak(1);
// 2.
exp(c)+exp(k)-(1-DELTAK)*exp(k(-1))+
PSI*(exp(k)/exp(k(-1))-1)^2=exp(a)*exp(k(-1))^(ALPHA);
// 3.
exp(i) = exp(k)-(1-DELTAK)*exp(k(-1));
// 4.
lambdak = 0;
// 5.
a = RHOA*a(-1)+erra;
kprev=k(-1);
end;
shocks;
var erra; stderr 0.015;
end;
steady;
stoch_simul(order=1,nocorr,nomoments,irf=0,print);