From c782ca368602fc1d57f83e19a82c5d53dfd3d76f Mon Sep 17 00:00:00 2001 From: Michel Juillard Date: Sun, 31 May 2015 16:48:20 +0200 Subject: [PATCH] adding occbin original files from 20140630 --- matlab/occbin/call_solve_one_constraint.m | 21 ++ matlab/occbin/call_solve_two_constraints.m | 21 ++ matlab/occbin/get_deriv.m | 83 +++++ matlab/occbin/get_pq.m | 28 ++ matlab/occbin/makechart.m | 81 +++++ matlab/occbin/makechart9.m | 136 ++++++++ matlab/occbin/map_regime.m | 25 ++ matlab/occbin/mkdata.m | 63 ++++ matlab/occbin/mkdatap_anticipated.m | 125 +++++++ .../occbin/mkdatap_anticipated_2constraints.m | 180 +++++++++++ matlab/occbin/pickaxes.m | 16 + matlab/occbin/process_constraint.m | 38 +++ matlab/occbin/setss.m | 14 + matlab/occbin/solve_no_constraint.m | 50 +++ matlab/occbin/solve_no_constraint_noclear.m | 48 +++ matlab/occbin/solve_one_constraint.1.m | 192 +++++++++++ matlab/occbin/solve_one_constraint.m | 200 ++++++++++++ matlab/occbin/solve_two_constraints.m | 305 ++++++++++++++++++ matlab/occbin/strmerge.m | 9 + matlab/occbin/tokenize.m | 55 ++++ tests/occbin/dynrbc.mod | 55 ++++ 21 files changed, 1745 insertions(+) create mode 100755 matlab/occbin/call_solve_one_constraint.m create mode 100755 matlab/occbin/call_solve_two_constraints.m create mode 100755 matlab/occbin/get_deriv.m create mode 100755 matlab/occbin/get_pq.m create mode 100755 matlab/occbin/makechart.m create mode 100755 matlab/occbin/makechart9.m create mode 100755 matlab/occbin/map_regime.m create mode 100755 matlab/occbin/mkdata.m create mode 100755 matlab/occbin/mkdatap_anticipated.m create mode 100755 matlab/occbin/mkdatap_anticipated_2constraints.m create mode 100755 matlab/occbin/pickaxes.m create mode 100755 matlab/occbin/process_constraint.m create mode 100755 matlab/occbin/setss.m create mode 100755 matlab/occbin/solve_no_constraint.m create mode 100755 matlab/occbin/solve_no_constraint_noclear.m create mode 100755 matlab/occbin/solve_one_constraint.1.m create mode 100755 matlab/occbin/solve_one_constraint.m create mode 100755 matlab/occbin/solve_two_constraints.m create mode 100755 matlab/occbin/strmerge.m create mode 100755 matlab/occbin/tokenize.m create mode 100755 tests/occbin/dynrbc.mod diff --git a/matlab/occbin/call_solve_one_constraint.m b/matlab/occbin/call_solve_one_constraint.m new file mode 100755 index 000000000..b4ec1baed --- /dev/null +++ b/matlab/occbin/call_solve_one_constraint.m @@ -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 diff --git a/matlab/occbin/call_solve_two_constraints.m b/matlab/occbin/call_solve_two_constraints.m new file mode 100755 index 000000000..b782e6db7 --- /dev/null +++ b/matlab/occbin/call_solve_two_constraints.m @@ -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 diff --git a/matlab/occbin/get_deriv.m b/matlab/occbin/get_deriv.m new file mode 100755 index 000000000..24f1b42ae --- /dev/null +++ b/matlab/occbin/get_deriv.m @@ -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 + diff --git a/matlab/occbin/get_pq.m b/matlab/occbin/get_pq.m new file mode 100755 index 000000000..5caab79af --- /dev/null +++ b/matlab/occbin/get_pq.m @@ -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; \ No newline at end of file diff --git a/matlab/occbin/makechart.m b/matlab/occbin/makechart.m new file mode 100755 index 000000000..a3801e956 --- /dev/null +++ b/matlab/occbin/makechart.m @@ -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 diff --git a/matlab/occbin/makechart9.m b/matlab/occbin/makechart9.m new file mode 100755 index 000000000..c6e65afcd --- /dev/null +++ b/matlab/occbin/makechart9.m @@ -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 diff --git a/matlab/occbin/map_regime.m b/matlab/occbin/map_regime.m new file mode 100755 index 000000000..a3ed88b47 --- /dev/null +++ b/matlab/occbin/map_regime.m @@ -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 + diff --git a/matlab/occbin/mkdata.m b/matlab/occbin/mkdata.m new file mode 100755 index 000000000..a7f855743 --- /dev/null +++ b/matlab/occbin/mkdata.m @@ -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); \ No newline at end of file diff --git a/matlab/occbin/mkdatap_anticipated.m b/matlab/occbin/mkdatap_anticipated.m new file mode 100755 index 000000000..524d3dd5d --- /dev/null +++ b/matlab/occbin/mkdatap_anticipated.m @@ -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,:); \ No newline at end of file diff --git a/matlab/occbin/mkdatap_anticipated_2constraints.m b/matlab/occbin/mkdatap_anticipated_2constraints.m new file mode 100755 index 000000000..b05c90d6d --- /dev/null +++ b/matlab/occbin/mkdatap_anticipated_2constraints.m @@ -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,:); \ No newline at end of file diff --git a/matlab/occbin/pickaxes.m b/matlab/occbin/pickaxes.m new file mode 100755 index 000000000..46c4b79c6 --- /dev/null +++ b/matlab/occbin/pickaxes.m @@ -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); diff --git a/matlab/occbin/process_constraint.m b/matlab/occbin/process_constraint.m new file mode 100755 index 000000000..d20911175 --- /dev/null +++ b/matlab/occbin/process_constraint.m @@ -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); \ No newline at end of file diff --git a/matlab/occbin/setss.m b/matlab/occbin/setss.m new file mode 100755 index 000000000..d1c0e440a --- /dev/null +++ b/matlab/occbin/setss.m @@ -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 + diff --git a/matlab/occbin/solve_no_constraint.m b/matlab/occbin/solve_no_constraint.m new file mode 100755 index 000000000..d9671409d --- /dev/null +++ b/matlab/occbin/solve_no_constraint.m @@ -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); + diff --git a/matlab/occbin/solve_no_constraint_noclear.m b/matlab/occbin/solve_no_constraint_noclear.m new file mode 100755 index 000000000..309c5c546 --- /dev/null +++ b/matlab/occbin/solve_no_constraint_noclear.m @@ -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); + diff --git a/matlab/occbin/solve_one_constraint.1.m b/matlab/occbin/solve_one_constraint.1.m new file mode 100755 index 000000000..34e296496 --- /dev/null +++ b/matlab/occbin/solve_one_constraint.1.m @@ -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_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 diff --git a/matlab/occbin/solve_one_constraint.m b/matlab/occbin/solve_one_constraint.m new file mode 100755 index 000000000..c446f4c40 --- /dev/null +++ b/matlab/occbin/solve_one_constraint.m @@ -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_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 diff --git a/matlab/occbin/solve_two_constraints.m b/matlab/occbin/solve_two_constraints.m new file mode 100755 index 000000000..dfca403a2 --- /dev/null +++ b/matlab/occbin/solve_two_constraints.m @@ -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_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_); + diff --git a/matlab/occbin/strmerge.m b/matlab/occbin/strmerge.m new file mode 100755 index 000000000..b8585df04 --- /dev/null +++ b/matlab/occbin/strmerge.m @@ -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 \ No newline at end of file diff --git a/matlab/occbin/tokenize.m b/matlab/occbin/tokenize.m new file mode 100755 index 000000000..f6e9f330c --- /dev/null +++ b/matlab/occbin/tokenize.m @@ -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 + diff --git a/tests/occbin/dynrbc.mod b/tests/occbin/dynrbc.mod new file mode 100755 index 000000000..51e612524 --- /dev/null +++ b/tests/occbin/dynrbc.mod @@ -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);