stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge branch 'master' into use-dynSeries

time-shift
Stéphane Adjemian (Charybdis) 2014-06-02 16:12:05 +02:00
parent 45d7f15c94 ff268520cd
commit a9e8dbb752
49 changed files with 2615 additions and 481 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

989
doc/dseries-and-reporting/dseriesReporting.tex
View File

File diff suppressed because it is too large Load Diff

125
doc/dynare.texi
View File

@ -10965,6 +10965,7 @@ ans =
@sp 1
@anchor{tex_rename}
@deftypefn{dseries} {@var{B} =} tex_rename (@var{A}, @var{name}, @var{newtexname})
Redefines the tex name of variable @var{name} to @var{newtexname}
@ -11037,13 +11038,15 @@ Computes yearly differences or growth rates.
@node Reporting
@chapter Reporting
Dynare provides a simple interface for creating @LaTeX{} reports,
comprised of @LaTeX{} tables and @code{PGFPLOTS/TikZ/PGF} graphs. You
can use the report as created through Dynare or pick out the pieces
(tables and graphs) you want for inclusion in your own paper. Though
Dynare provides a subset of options available through @code{PGFPLOTS},
you can easily modify the graphs created by Dynare using the options
available in the @code{PGFPLOTS} manual.
Dynare provides a simple interface for creating @LaTeX{} reports, comprised of
@LaTeX{} tables and @code{PGFPLOTS/Ti}@i{k}@code{Z} graphs. You can use the
report as created through Dynare or pick out the pieces (tables and graphs) you
want for inclusion in your own paper. Though Dynare provides a subset of
options available through @code{PGFPLOTS/Ti}@i{k}@code{Z}, you can easily
modify the graphs created by Dynare using the options available in the
@code{PGFPLOTS/Ti}@i{k}@code{Z} manual. You can either do this manually or by
passing the options to @ref{miscTikzAxisOptions}, @ref{miscTikzAxisOptions}, or
@ref{graphMiscTikzAddPlotOptions}.
Reports are created and modified by calling methods on class
objects. The objects are hierarchical, with the following order (from
@ -11159,7 +11162,8 @@ command. Default: @code{`!'}
@end table
@end defmethod
@defmethod Report addGraph data, graphDirName, graphName, graphSize, height, showGrid, showLegend, showLegendBox, legendLocation, legendOrientation, legendFontSize, seriesToUse, shade, shadeColor, shadeOpacity, title, width, xlabel, ylabel, xAxisTight, xrange, xTicks, xTickLabels, xTickLabelAnchor, xTickLabelRotation, yAxisTight, yrange, showZeroline
@anchor{addGraph}
@defmethod Report addGraph data, graphDirName, graphName, graphSize, height, showGrid, showLegend, showLegendBox, legendLocation, legendOrientation, legendFontSize, miscTikzAxisOptions, miscTikzPictureOptions, seriesToUse, shade, shadeColor, shadeOpacity, title, width, xlabel, ylabel, xAxisTight, xrange, xTicks, xTickLabels, xTickLabelAnchor, xTickLabelRotation, yAxisTight, yrange, showZeroline
Adds a @code{Graph} to a @code{Section}.
@optionshead
@table @code
@ -11189,8 +11193,12 @@ The height of the graph, in inches. Default: @code{4.5}
Whether or not to display the major grid on the graph. Default:
@code{true}
@anchor{showLegend}
@item showLegend, @code{BOOLEAN}
Whether or not to display the legend. Default: @code{false}
Whether or not to display the legend. NB: Unless you use the
@ref{graphLegendName} option, the name displayed in the legend is the
@code{tex} name associated with the @code{dseries}. You can modify this
@code{tex} name by using @ref{tex_rename}. Default: @code{false}
@item showLegendBox, @code{BOOLEAN}
Whether or not to display a box around the legend. Default:
@ -11205,6 +11213,24 @@ Orientation of the legend. Default: @code{`horizontal'}
@item legendFontSize, @code{`tiny'} | @code{`scriptsize'} | @code{`footnotesize'} | @code{`small'} | @code{`normalsize'} | @code{`large'} | @code{`Large'} | @code{`LARGE'} | @code{`huge'} | @code{`Huge'}
The font size for legend entries. Default: @code{tiny}
@anchor{miscTikzAxisOptions}
@item miscTikzAxisOptions, @code{STRING}
If you are comfortable with @code{PGFPLOTS/Ti}@i{k}@code{Z}, you can use this
option to pass arguments directly to the @code{PGFPLOTS/Ti}@i{k}@code{Z}
@code{axis} environment command. Specifically to be used for desired
@code{PGFPLOTS/Ti}@i{k}@code{Z} options that have not been incorporated into
Dynare Reproting. Default: @code{empty}
@anchor{miscTikzPictureOptions}
@item miscTikzPictureOptions, @code{STRING}
If you are comfortable with @code{PGFPLOTS/Ti}@i{k}@code{Z}, you can use this
option to pass arguments directly to the @code{PGFPLOTS/Ti}@i{k}@code{Z}
@code{tikzpicture} environment command. (@i{e.g.,} to scale the graph in the x
and y dimensions, you can pass following to this option: @code{`xscale=2.5,
yscale=0.5'}). Specifically to be used for desired
@code{PGFPLOTS/Ti}@i{k}@code{Z} options that have not been incorporated into
Dynare Reproting. Default: @code{empty}
@anchor{seriesToUse}
@item seriesToUse, @code{CELL_ARRAY_STRINGS}
The names of the series contained in the @code{dseries} provided to
@ -11242,8 +11268,8 @@ The x-axis label. Default: @code{none}
The y-axis label. Default: @code{none}
@item xAxisTight, @code{BOOLEAN}
Use a tight x axis. If false, uses pgfplots @code{enlarge x limits} to
choose appropriate axis size. Default: @code{true}
Use a tight x axis. If false, uses @code{PGFPLOTS/Ti}@i{k}@code{Z}
@code{enlarge x limits} to choose appropriate axis size. Default: @code{true}
@item xrange, @code{dates}
The boundary on the x-axis to display in the graph. Default: all
@ -11269,8 +11295,8 @@ Where to anchor the x tick label. Default: @code{`south'}
The amount to rotate the x tick labels by. Default: @code{0}
@item yAxisTight, @code{BOOLEAN}
Use a tight y axis. If false, uses pgfplots @code{enlarge y limits} to
choose appropriate axis size. Default: @code{false}
Use a tight y axis. If false, uses @code{PGFPLOTS/Ti}@i{k}@code{Z}
@code{enlarge y limits} to choose appropriate axis size. Default: @code{false}
@item yrange, @code{NUMERICAL_VECTOR}
The boundary on the y-axis to display in the graph, represented as a
@ -11331,14 +11357,25 @@ Whether or not to show vertical lines separating the columns. Default: @code{fal
@end defmethod
@anchor{addSeries}
@defmethod Report addSeries data, graphLineColor, graphLineStyle, graphLineWidth, graphMarker, graphMarkerEdgeColor, graphMarkerFaceColor, graphMarkerSize, tableDataRhs, tableRowColor, tableRowIndent, tableShowMarkers, tableAlignRight, tableNegColor, tablePosColor, tableSubSectionHeader, zeroTol
Adds a @code{Series} to a @code{Graph} or a @code{Table}.
@defmethod Report addSeries data, graphHline, graphLegendName, graphLineColor, graphLineStyle, graphLineWidth, graphMarker, graphMarkerEdgeColor, graphMarkerFaceColor, graphMarkerSize, graphMiscTikzAddPlotOptions, graphShowInLegend, graphVline, tableDataRhs, tableRowColor, tableRowIndent, tableShowMarkers, tableAlignRight, tableNegColor, tablePosColor, tableSubSectionHeader, zeroTol
Adds a @code{Series} to a @code{Graph} or a @code{Table}. NB: Options specific
to graphs begin with `@code{graph}' while options specific to tables begin with
`@code{table}'.
@optionshead
@table @code
@item data, @code{dseries}
@xref{data}.
@item graphHline, @code{DOUBLE}
Use this option to draw a horizontal line at the given value. Default:
@code{empty}
@anchor{graphLegendName}
@item graphLegendName, @code{STRING}
The name to display in the legend for this series. Will be displayed only if
the @ref{data} option has been set. Default: the @code{tex} name of the series
@item graphLineColor, @code{`red'} | @code{`green'} | @code{`blue'} | @code{`cyan'} | @code{`magenta'} | @code{`yellow'} | @code{`black'} | @code{`gray'} | @code{`darkgray'} | @code{`lightgray'} | @code{`brown'} | @code{`lime'} | @code{`olive'} | @code{`orange'} | @code{`pink'} | @code{`purple'} | @code{`teal'} | @code{`violet'} | @code{`white'}
Color to use for the series in a graph. Default: @code{`black'}
@ -11360,6 +11397,23 @@ The face color of the graph marker. Default: @code{graphLineColor}
@item graphMarkerSize, @code{DOUBLE}
The size of the graph marker. Default: @code{1}
@anchor{graphMiscTikzAddPlotOptions}
@item graphMiscTikzAddPlotOptions, @code{STRING}
If you are comfortable with @code{PGFPLOTS/Ti}@i{k}@code{Z}, you can use this
option to pass arguments directly to the @code{PGFPLOTS/Ti}@i{k}@code{Z}
@code{addPlots} command. (@i{e.g.,} Instead of passing the marker options
above, you can pass a string such as the following to this option:
@code{`mark=halfcircle*,mark options=@{rotate=90,scale=3@}'}). Specifically to be
used for desired @code{PGFPLOTS/Ti}@i{k}@code{Z} options that have not been
incorporated into Dynare Reproting. Default: @code{empty}
@item graphShowInLegend, @code{BOOLEAN}
Whether or not to show this series in the legend, given that the
@ref{showLegend} option was passed to @ref{addGraph}. Default: @code{true}
@item graphVline, @code{dates}
Use this option to draw a vertical line at a given date. Default: @code{empty}
@item tableDataRhs, @code{dseries}
A series to be added to the right of the current series. Usefull for
displaying aggregate data for a series. @i{e.g} if the series is
@ -11406,9 +11460,36 @@ with it. It is equivalent to adding an empty series with a
name. Default: @code{''}
@item zeroTol, @code{DOUBLE}
The zero tolerance. Anything smaller than @code{zeroTol} and larger
than @code{-zeroTol} will be set to zero before being
graphed. Default: @math{1e-6}
The zero tolerance. Anything smaller than @code{zeroTol} and larger than
@code{-zeroTol} will be set to zero before being graphed or written to the
table. Default: @math{1e-6}
@end table
@end defmethod
@defmethod Report addParagraph balancedCols, cols, heading, indent, text
Adds a @code{Paragraph} to a @code{Section}. NB: The @code{Section} can only be
comprised of @code{Paragraphs} and must only have 1 column.
@optionshead
@table @code
@item balancedCols, @code{BOOLEAN}
Determines whether the text is spread out evenly across the columns when the
@code{Paragraph} has more than one column. Default: @code{true}
@item cols, @code{INTEGER}
The number of columns for the @code{Paragraph}. Default: @code{1}
@item heading, @code{STRING}
The heading for the Paragraph (like a section heading). The string must be
correct @LaTeX{} code. Default: @code{empty}
@item indent, @code{BOOLEAN}
Whether or not to indent the paragraph. Default: @code{true}
@item text, @code{STRING}
The paragraph itself. The string must be correct @LaTeX{} code. Default:
@code{empty}
@end table
@end defmethod
@ -11468,16 +11549,16 @@ rep = rep.addSection(`cols', 2);
rep = rep.addGraph(`title', `Graph (1,1)', `showLegend', true, ...
`xrange', dates(`2007q1'):dates(`2013q4'), ...
`shade', dates(`2012q2'):dates(`2013q4'));
rep = rep.addSeries(`data', dsq@{`SERIES1'@}, `color', `b', ...
rep = rep.addSeries(`data', dsq@{`SERIES1'@}, `graphLineColor', `blue', ...
`graphLineWidth', 1);
rep = rep.addSeries(`data', dsq@{`SERIES2'@}, `color', `g', ...
rep = rep.addSeries(`data', dsq@{`SERIES2'@}, `graphLineColor', `green', ...
`graphLineStyle', '--', `graphLineWidth', 1.5);
rep = rep.addGraph(`title', `Graph (1,2)', `showLegend', true, ...
`xrange', dates(`2007q1'):dates(`2013q4'), ...
`shade', dates(`2012q2'):dates(`2013q4'));
rep = rep.addSeries(`data', dsq@{`SERIES3'@}, `color', `b', ...
rep = rep.addSeries(`data', dsq@{`SERIES3'@}, `graphLineColor', `blue', ...
`graphLineWidth', 1);
rep = rep.addSeries(`data', dsq@{`SERIES4'@}, `color', `g', ...
rep = rep.addSeries(`data', dsq@{`SERIES4'@}, `graphLineColor', `green', ...
`graphLineStyle', '--', `graphLineWidth', 1.5);
% Section 2

1
matlab/dynare_config.m
View File

@ -57,6 +57,7 @@ addpath([dynareroot '/parallel/'])
addpath([dynareroot '/particle/'])
addpath([dynareroot '/gsa/'])
addpath([dynareroot '/ep/'])
addpath([dynareroot '/lmmcp/'])
addpath([dynareroot '/utilities/doc/'])
addpath([dynareroot '/utilities/tests/'])
addpath([dynareroot '/utilities/dates/'])

8
matlab/dynare_solve.m
View File

@ -164,6 +164,14 @@ elseif options_.solve_algo == 3
else
[x,info] = csolve(func,x,[],1e-6,500,varargin{:});
end
elseif options_.solve_algo == 10
olmmcp = options_.lmmcp;
[x,fval,exitflag] = lmmcp(func,x,olmmcp.lb,olmmcp.ub,olmmcp,varargin{:});
if exitflag == 1
info = 0;
else
info = 1;
end
else
error('DYNARE_SOLVE: option solve_algo must be one of [0,1,2,3,4,9]')
end

33
matlab/ep/ep_problem_2.m
View File

@ -26,6 +26,7 @@ i_cols_1 = pm.i_cols_1;
i_cols_j = pm.i_cols_j;
icA = pm.icA;
i_cols_T = pm.i_cols_T;
eq_index = pm.eq_index;
i_cols_p = i_cols(1:nyp);
i_cols_s = i_cols(nyp+(1:ny));
@ -92,15 +93,18 @@ for i = 1:order+1
% in first period we don't keep track of
% predetermined variables
i_cols_A = [i_cols_As - ny; i_cols_Af];
A1(i_rows,i_cols_A) = A1(i_rows,i_cols_A) + weights(k)*jacobian(:,i_cols_1);
A1(i_rows,i_cols_A) = A1(i_rows,i_cols_A) + weights(k)*jacobian(eq_index,i_cols_1);
else
i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af];
A1(i_rows,i_cols_A) = A1(i_rows,i_cols_A) + weights(k)*jacobian(:,i_cols_j);
A1(i_rows,i_cols_A) = A1(i_rows,i_cols_A) + weights(k)*jacobian(eq_index,i_cols_j);
end
else
d1 = dynamic_model(z,innovation,params,steady_state,i+1);
end
res(:,i,1) = res(:,i,1)+weights(k)*d1;
if any(isnan(d1))
pause
end
res(:,i,1) = res(:,i,1)+weights(k)*d1(eq_index);
end
if nargout > 1
[ir,ic,v] = find(A1);
@ -121,12 +125,15 @@ for i = 1:order+1
if nargout > 1
[d1,jacobian] = dynamic_model(z,innovation,params,steady_state,i+1);
i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af];
[ir,ic,v] = find(jacobian(:,i_cols_j));
[ir,ic,v] = find(jacobian(eq_index,i_cols_j));
nzA{i,j} = [i_rows(ir),i_cols_A(ic), v]';
else
d1 = dynamic_model(z,innovation,params,steady_state,i+1);
end
res(:,i,j) = d1;
if any(isnan(d1))
pause
end
res(:,i,j) = d1(eq_index);
if nargout > 1
i_cols_Af = i_cols_Af + ny;
end
@ -142,13 +149,16 @@ for i = 1:order+1
if nargout > 1
[d1,jacobian] = dynamic_model(z,innovation,params,steady_state,i+1);
i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af];
[ir,ic,v] = find(jacobian(:,i_cols_j));
[ir,ic,v] = find(jacobian(eq_index,i_cols_j));
nzA{i,j} = [i_rows(ir),i_cols_A(ic),v]';
i_cols_Af = i_cols_Af + ny;
else
d1 = dynamic_model(z,innovation,params,steady_state,i+1);
end
res(:,i,j) = d1;
if any(isnan(d1))
pause
end
res(:,i,j) = d1(eq_index);
end
i_rows = i_rows + ny;
if mod(j,nnodes) == 0
@ -172,16 +182,19 @@ for j=1:world_nbr
[d1,jacobian] = dynamic_model(Y(i_rows_y,j),x,params, ...
steady_state,i+1);
if i < periods
[ir,ic,v] = find(jacobian(:,i_cols_j));
[ir,ic,v] = find(jacobian(eq_index,i_cols_j));
else
[ir,ic,v] = find(jacobian(:,i_cols_T));
[ir,ic,v] = find(jacobian(eq_index,i_cols_T));
end
nzA{i,j} = [offset_r+ir,offset_c+icA(ic), v]';
else
d1 = dynamic_model(Y(i_rows_y,j),x,params, ...
steady_state,i+1);
end
res(:,i,j) = d1;
if any(isnan(d1))
pause
end
res(:,i,j) = d1(eq_index);
i_rows_y = i_rows_y + ny;
offset_c = offset_c + world_nbr*ny;
offset_r = offset_r + world_nbr*ny;

8
matlab/ep/extended_path.m
View File

@ -150,6 +150,14 @@ if options_.ep.stochastic.order > 0
pfm.nnodes = nnodes;
end
% compute number of blocks
[block_nbr,pfm.world_nbr] = get_block_world_nbr(options_.ep.stochastic.algo,nnodes,options_.ep.ut.k,options_.ep.periods);
% set boundaries if mcp
[lb,ub,pfm.eq_index] = get_complementarity_conditions(M_);
options_.lmmcp.lb = repmat(lb,block_nbr,1);
options_.lmmcp.ub = repmat(ub,block_nbr,1);
pfm.block_nbr = block_nbr;
% Main loop.
while (t<sample_size)
if ~mod(t,10)

31
matlab/ep/get_block_world_nbr.m Normal file
View File

@ -0,0 +1,31 @@
function [block_nbr,world_nbr] = get_block_world_nbr(algo,nnodes,order,periods)
% Copyright (C) 2014 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
switch algo
case 0
world_nbr = nnodes^order;
block_nbr = 1+(nnodes^(order+1)-nnodes)/(nnodes-1)+(periods-order)*world_nbr;
case 1
world_nbr = 1+(nnodes-1)*order;
block_nbr = (order+(nnodes-1)*(order-1)*order/2+(periods-order)* ...
world_nbr);
otherwise
error('This case is not supposed to happen')
end

11
matlab/ep/solve_stochastic_perfect_foresight_model_1.m
View File

@ -84,7 +84,7 @@ end
% The second row block is ny x nnodes
% The third row block is ny x nnodes^2
% and so on until size ny x nnodes^order
world_nbr = 1+(nnodes-1)*order;
world_nbr = pfm.world_nbr;
Y = endo_simul(:,2:end-1);
Y = repmat(Y,1,world_nbr);
pfm.y0 = endo_simul(:,1);
@ -92,10 +92,9 @@ pfm.y0 = endo_simul(:,1);
% The columns of A map the elements of Y such that
% each block of Y with ny rows are unfolded column wise
% number of blocks
block_nbr = (order+(nnodes-1)*(order-1)*order/2+(periods-order)*world_nbr);
block_nbr = pfm.block_nbr;
% dimension of the problem
dimension = ny*block_nbr;
pfm.block_nbr = block_nbr;
pfm.dimension = dimension;
if order == 0
i_upd_r = (1:ny*periods);
@ -143,11 +142,13 @@ pfm.i_cols_T = i_cols_T;
pfm.i_upd_r = i_upd_r;
pfm.i_upd_y = i_upd_y;
options_.solve_algo = 9;
options_.steady.maxit = 100;
y = repmat(steady_state,block_nbr,1);
old_options = options_;
options_.solve_algo = options_.ep.solve_algo;
options_.steady.maxit = options_.ep.maxit;
[y,info] = dynare_solve(@ep_problem_2,y,1,exo_simul,pfm);
options_ = old_options;
if info
flag = 1;
err = info;

148
matlab/lmmcp/catstruct.m Normal file
View File

@ -0,0 +1,148 @@
function A = catstruct(varargin)
% CATSTRUCT Concatenate or merge structures with different fieldnames
% X = CATSTRUCT(S1,S2,S3,...) merges the structures S1, S2, S3 ...
% into one new structure X. X contains all fields present in the various
% structures. An example:
%
% A.name = 'Me' ;
% B.income = 99999 ;
% X = catstruct(A,B)
% % -> X.name = 'Me' ;
% % X.income = 99999 ;
%
% If a fieldname is not unique among structures (i.e., a fieldname is
% present in more than one structure), only the value from the last
% structure with this field is used. In this case, the fields are
% alphabetically sorted. A warning is issued as well. An axample:
%
% S1.name = 'Me' ;
% S2.age = 20 ; S3.age = 30 ; S4.age = 40 ;
% S5.honest = false ;
% Y = catstruct(S1,S2,S3,S4,S5) % use value from S4
%
% The inputs can be array of structures. All structures should have the
% same size. An example:
%
% C(1).bb = 1 ; C(2).bb = 2 ;
% D(1).aa = 3 ; D(2).aa = 4 ;
% CD = catstruct(C,D) % CD is a 1x2 structure array with fields bb and aa
%
% The last input can be the string 'sorted'. In this case,
% CATSTRUCT(S1,S2, ..., 'sorted') will sort the fieldnames alphabetically.
% To sort the fieldnames of a structure A, you could use
% CATSTRUCT(A,'sorted') but I recommend ORDERFIELDS for doing that.
%
% When there is nothing to concatenate, the result will be an empty
% struct (0x0 struct array with no fields).
%
% NOTE: To concatenate similar arrays of structs, you can use simple
% concatenation:
% A = dir('*.mat') ; B = dir('*.m') ; C = [A ; B] ;
%
% See also CAT, STRUCT, FIELDNAMES, STRUCT2CELL, ORDERFIELDS
% for Matlab R13 and up
% version 3.0 (mar 2013)
% (c) Jos van der Geest
% email: jos@jasen.nl
% (C) 2013 Christophe Gouel
% History
% Created in 2005
% Revisions
% 2.0 (sep 2007) removed bug when dealing with fields containing cell
% arrays (Thanks to Rene Willemink)
% 2.1 (sep 2008) added warning and error identifiers
% 2.2 (oct 2008) fixed error when dealing with empty structs (Thanks to
% Lars Barring)
% 3.0 (mar 2013) fixed problem when the inputs were array of structures
% (thanks to Tor Inge Birkenes for pointing this out).
% Rephrased the help section as well.
error(nargchk(1,Inf,nargin)) ;
N = nargin ;
if ~isstruct(varargin{end}),
if isequal(varargin{end},'sorted'),
sorted = 1 ;
N = N-1 ;
error(nargchk(1,Inf,N)) ;
else
error('catstruct:InvalidArgument','Last argument should be a structure, or the string "sorted".') ;
end
else
sorted = 0 ;
end
sz0 = [] ; % used to check that all inputs have the same size
% used to check for a few trivial cases
NonEmptyInputs = false(N,1) ;
NonEmptyInputsN = 0 ;
% used to collect the fieldnames and the inputs
FN = cell(N,1) ;
VAL = cell(N,1) ;
% parse the inputs
for ii=1:N,
X = varargin{ii} ;
if ~isstruct(X),
error('catstruct:InvalidArgument',['Argument #' num2str(ii) ' is not a structure.']) ;
end
if ~isempty(X),
% empty structs are ignored
if ii > 1 && ~isempty(sz0)
if ~isequal(size(X), sz0)
error('catstruct:UnequalSizes','All structures should have the same size.') ;
end
else
sz0 = size(X) ;
end
NonEmptyInputsN = NonEmptyInputsN + 1 ;
NonEmptyInputs(ii) = true ;
FN{ii} = fieldnames(X) ;
VAL{ii} = struct2cell(X) ;
end
end
if NonEmptyInputsN == 0
% all structures were empty
A = struct([]) ;
elseif NonEmptyInputsN == 1,
% there was only one non-empty structure
A = varargin{NonEmptyInputs} ;
if sorted,
A = orderfields(A) ;
end
else
% there is actually something to concatenate
FN = cat(1,FN{:}) ;
VAL = cat(1,VAL{:}) ;
FN = squeeze(FN) ;
VAL = squeeze(VAL) ;
MatlabVersion = version;
if isoctave || str2double(MatlabVersion(end-5:end-2))<2013 % Equivalent to, but faster than if verLessThan('matlab','8.1')
[UFN,ind] = unique(FN) ;
else
[UFN,ind] = unique(FN,'legacy') ;
end
if numel(UFN) ~= numel(FN),
warning('catstruct:DuplicatesFound','Fieldnames are not unique between structures.') ;
sorted = 1 ;
end
if sorted,
VAL = VAL(ind,:) ;
FN = FN(ind,:) ;
end
A = cell2struct(VAL, FN);
A = reshape(A, sz0) ; % reshape into original format
end

77
matlab/lmmcp/dyn_lmmcp.m Normal file
View File

@ -0,0 +1,77 @@
function [endo_simul,info] = dyn_lmmcp(M,options,oo)
% Copyright (C) 2014 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
[lb,ub,eq_index] = get_complementarity_conditions(M);
lead_lag_incidence = M.lead_lag_incidence;
ny = M.endo_nbr;
max_lag = M.maximum_endo_lag;
nyp = nnz(lead_lag_incidence(1,:)) ;
iyp = find(lead_lag_incidence(1,:)>0) ;
ny0 = nnz(lead_lag_incidence(2,:)) ;
iy0 = find(lead_lag_incidence(2,:)>0) ;
nyf = nnz(lead_lag_incidence(3,:)) ;
iyf = find(lead_lag_incidence(3,:)>0) ;
nd = nyp+ny0+nyf;
nrc = nyf+1 ;
isp = [1:nyp] ;
is = [nyp+1:ny+nyp] ;
isf = iyf+nyp ;
isf1 = [nyp+ny+1:nyf+nyp+ny+1] ;
stop = 0 ;
iz = [1:ny+nyp+nyf];
periods = options.periods;
steady_state = oo.steady_state;
params = M.params;
endo_simul = oo.endo_simul;
exo_simul = oo.exo_simul;
i_cols_1 = nonzeros(lead_lag_incidence(2:3,:)');
i_cols_A1 = find(lead_lag_incidence(2:3,:)');
i_cols_T = nonzeros(lead_lag_incidence(1:2,:)');
i_cols_j = 1:nd;
i_upd = ny+(1:periods*ny);
x = endo_simul(:);
model_dynamic = str2func([M.fname,'_dynamic']);
z = x(find(lead_lag_incidence'));
[res,A] = model_dynamic(z, exo_simul, params, steady_state,2);
nnzA = nnz(A);
LB = repmat(lb,periods,1);
UB = repmat(ub,periods,1);
Y0 = endo_simul(:,1);
YT = endo_simul(:,end);
x = endo_simul(:,2:end-1);
x = x(:);
func_handle = @(x) dyn_lmmcp_func(x,model_dynamic, Y0, YT, exo_simul, ...
params, steady_state, periods, ny, ...
lead_lag_incidence, i_cols_A1, i_cols_1, ...
i_cols_T, i_cols_j,nnzA,eq_index);
[x, info] = lmmcp(func_handle,x,LB,UB);
endo_simul = [Y0 reshape(x,ny,periods) YT];

56
matlab/lmmcp/dyn_lmmcp_func.m Normal file
View File

@ -0,0 +1,56 @@
function [F,A] = dyn_lmmcp_func(x, model_dynamic, Y0, YT, exo_simul, params, ...
steady_state, periods, ny, lead_lag_incidence, ...
i_cols_A1, i_cols_1, i_cols_T, i_cols_j, ...
nnzA,eq_index)
% Copyright (C) 2014 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
Y = [Y0; x; YT];
F = zeros(periods*ny,1);
if nargout == 2
A = sparse([],[],[],periods*ny,periods*ny,periods*nnzA);
end
i_rows = 1:ny;
i_cols = find(lead_lag_incidence');
i_cols_A = i_cols;
for it = 2:(periods+1)
[res,jacobian] = model_dynamic(Y(i_cols),exo_simul, params, ...
steady_state,it);
F(i_rows) = res(eq_index);
if nargout == 2
if it == 2
A(i_rows,i_cols_A1) = jacobian(eq_index,i_cols_1);
elseif it == periods+1
A(i_rows,i_cols_A(i_cols_T)) = jacobian(eq_index,i_cols_T);
else
A(i_rows,i_cols_A) = jacobian(eq_index,i_cols_j);
end
end
i_rows = i_rows + ny;
i_cols = i_cols + ny;
if nargout == 2 && it > 2
i_cols_A = i_cols_A + ny;
end
end

56
matlab/lmmcp/get_complementarity_conditions.m Normal file
View File

@ -0,0 +1,56 @@
function [lb,ub,eq_index] = get_complementarity_conditions(M)
% Copyright (C) 2014 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
etags = M.equations_tags;
ub = inf(M.endo_nbr,1);
lb = -ub;
eq_index = (1:M.endo_nbr)';
for i=1:size(etags,1)
if strcmp(etags{i,2},'mcp')
str = etags{i,3};
kop = strfind(etags{i,3},'<');
if ~isempty(kop)
k = find(strcmp(strtrim(str(1:kop-1)),cellstr(M.endo_names)));
if isempty(k)
error(sprintf(['Complementarity condition %s: variable %s is ' ...
'not recognized',etags{i,3},b{1}]))
end
ub(k) = str2num(str(kop+1:end));
eq_index(etags{i,1}) = k;
eq_index(k) = etags{i,1};
else
kop = strfind(etags{i,3},'>');
if ~isempty(kop)
k = find(strcmp(strtrim(str(1:kop-1)),cellstr(M.endo_names)));
if isempty(k)
error(sprintf(['Complementarity condition %s: variable %s is ' ...
'not recognized',etags{i},b{1}]))
end
lb(k) = str2num(str(kop+1:end));
eq_index(etags{i,1}) = k;
eq_index(k) = etags{i,1};
else
error(sprintf(['Complementarity condition %s can''t be ' ...
'parsed'],etags{i,3}))
end
end
end
end

614
matlab/lmmcp/lmmcp.m Normal file
View File

@ -0,0 +1,614 @@
function [x,FVAL,EXITFLAG,OUTPUT,JACOB] = lmmcp(FUN,x,lb,ub,options,varargin)
% LMMCP solves a mixed complementarity problem.
%
% LMMCP uses a semismooth least squares formulation. The method applies a
% Levenberg-Marquardt/Gauss-Newton algorithm to a least-squares formulation.
%
% X = LMMCP(FUN,X0) tries to solve the system of nonlinear equations F(X)=0 and
% starts at the vector X0. FUN accepts a vector X and return a vector F of equation
% values F evaluated at X and, as second output if required, a matrix J, the
% Jacobian evaluated at X.
%
% X = LMMCP(FUN,X0,LB,UB) solves the mixed complementarity problem of the form:
% LB =X => F(X)>0,
% LB<=X<=UB => F(X)=0,
% X =UB => F(X)<0.
%
% X = LMMCP(FUN,X0,LB,UB,OPTIONS) solves the MCP problem using the options
% defined in the structure OPTIONS. Main fields are
% Display : control the display of iterations, 'none' (default),
% 'iter-detailed' or 'final-detailed'
% Switch from phase I to phase II
% preprocess : activate preprocessor for phase I (default = 1)
% presteps : number of iterations in phase I (default = 20)
% Termination parameters
% MaxIter : Maximum number of iterations (default = 500)
% tmin : safeguard stepsize (default = 1E-12)
% TolFun : Termination tolerance on the function value, a positive
% scalar (default = sqrt(eps))
% Stepsize parameters
% m : number of previous function values to use in the nonmonotone
% line search rule (default = 10)
% kwatch : maximum number of steps (default = 20 and should not be
% smaller than m)
% watchdog : activate the watchdog strategy (default = 1)
% Ther are other minor parameters. Please see the code for their default values
% and interpretation.
%
% [X,FVAL] = LMMCP(FUN,X0,...) returns the value of the equations FUN at X.
%
% [X,FVAL,EXITFLAG] = LMMCP(FUN,X0,...) returns EXITFLAG that describes the exit
% conditions. Possible values are
% 1 : LMMCP converged to a root
% 0 : Too many iterations
% -1 :
%
% [X,FVAL,EXITFLAG,OUTPUT] = LMMCP(FUN,X0,...) returns the structure OUTPUT that
% contains the number of iterations (OUTPUT.iterations), the value of the merit
% function (OUTPUT.Psix), and the norm of the derivative of the merit function
% (OUTPUT.normDPsix).
%
% [X,FVAL,EXITFLAG,OUTPUT,JACOB] = LMMCP(FUN,X0,...) returns JACOB the Jacobian
% of FUN evaluated at X.
%
% More details of the main program may be found in the following paper:
%
% Christian Kanzow and Stefania Petra: On a semismooth least squares formulation of
% complementarity problems with gap reduction. Optimization Methods and Software
% 19, 2004, pp. 507-525.
%
% In addition, the current implementation uses a preprocessor which is the
% projected Levenberg-Marquardt step from the following preprint:
%
% Christian Kanzow and Stefania Petra: Projected filter trust region methods for a
% semismooth least squares formulation of mixed complementarity
% problems. Optimization Methods and Software
% 22, 2007, pp. 713-735.
%
% A user's guide is also available:
%
% Christian Kanzow and Stefania Petra (2005).
% LMMCP --- A Levenberg-Marquardt-type MATLAB Solver for Mixed Complementarity Problems.
% University of Wuerzburg.
% http://www.mathematik.uni-wuerzburg.de/~kanzow/software/UserGuide.pdf
%
% This is a modification by Christophe Gouel of the original files, which can be
% downloaded from:
% http://www.mathematik.uni-wuerzburg.de/~kanzow/software/LMMCP.zip
%
% copyright: Christian Kanzow and Stefania Petra
% Institute of Applied Mathematics and Statistics
% University of Wuerzburg
% Am Hubland
% 97074 Wuerzburg
% GERMANY
%
% e-mail: kanzow@mathematik.uni-wuerzburg.de
% petra@mathematik.uni-wuerzburg.de
%% Initialization
defaultopt = struct(...
'beta', 0.55,...
'Big', 1e10,...
'delta', 5,...
'deltamin', 1,...
'Display', 'none',...
'epsilon1', 1e-6,...
'eta', 0.95,...
'kwatch', 20,...
'lambda1', 0.1,...
'm', 10,...
'MaxIter', 500,...
'null', 1e-10,...
'preprocess', 1,...
'presteps', 20,...
'sigma', 1e-4,...
'sigma1', 0.5,...
'sigma2', 2,...
'tmin', 1e-12,...
'TolFun', sqrt(eps),...
'watchdog', 1);
if nargin < 4
ub = inf(size(x));
if nargin < 3
lb = -inf(size(x));
end
end
if nargin < 5 || isempty(options) || ~isstruct(options)
options = defaultopt;
else
warning('off','catstruct:DuplicatesFound')
options = catstruct(defaultopt,options);
end
warning('off','MATLAB:rankDeficientMatrix')
switch options.Display
case {'off','none'}
verbosity = 0;
case {'iter','iter-detailed'}
verbosity = 2;
case {'final','final-detailed'}
verbosity = 1;
otherwise
verbosity = 0;
end
% parameter settings
eps1 = options.epsilon1;
eps2 = 0.5*options.TolFun^2;
null = options.null;
Big = options.Big;
% maximal number of iterations
kmax = options.MaxIter;
% choice of lambda
lambda1 = options.lambda1;
lambda2 = 1-lambda1;
% steplength parameters
beta = options.beta;
sigma = options.sigma;
tmin = options.tmin;
% parameters watchdog and nonmonotone line search; redefined later
m = options.m;
kwatch = options.kwatch;
watchdog = options.watchdog; % 1=watchdog strategy active, otherwise not
% parameters for preprocessor
preprocess = options.preprocess; % 1=preprocessor used, otherwise not
presteps = options.presteps; % maximum number of preprocessing steps
% trust-region parameters for preprocessor
delta = options.delta;
deltamin = options.deltamin;
sigma1 = options.sigma1;
sigma2 = options.sigma2;
eta = options.eta;
% initializations
k = 0;
k_main = 0;
% compute a feasible starting point by projection
x = max(lb,min(x,ub));
n = length(x);
OUTPUT.Dim = n;
% definition of index sets I_l, I_u and I_lu
Indexset = zeros(n,1);
I_l = lb>-Big & ub>Big;
I_u = lb<-Big & ub<Big;
I_lu = lb>-Big & ub<Big;
Indexset(I_l) = 1;
Indexset(I_u) = 2;
Indexset(I_lu) = 3;
% function evaluations
[Fx,DFx] = feval(FUN,x,varargin{:});
% choice of NCP-function and corresponding evaluations
Phix = Phi(x,Fx,lb,ub,lambda1,lambda2,n,Indexset);
normPhix = norm(Phix);
Psix = 0.5*(Phix'*Phix);
DPhix = DPhi(x,Fx,DFx,lb,ub,lambda1,lambda2,n,Indexset);
DPsix = DPhix'*Phix;
normDPsix = norm(DPsix);
% save initial values
x0 = x;
Phix0 = Phix;
Psix0 = Psix;
DPhix0 = DPhix;
DPsix0 = DPsix;
normDPsix0 = normDPsix;
% watchdog strategy
aux = zeros(m,1);
aux(1) = Psix;
MaxPsi = Psix;
if watchdog==1
kbest = k;
xbest = x;
Phibest = Phix;
Psibest = Psix;
DPhibest = DPhix;
DPsibest = DPsix;
normDPsibest = normDPsix;
end
% initial output
if verbosity > 1
fprintf(' k Psi(x) || DPsi(x) || stepsize\n');
disp('====================================================================')
disp('********************* Output at starting point *********************')
fprintf('%4.0f %24.5e %24.5e\n',k,Psix,normDPsix);
end
%% Preprocessor using local method
if preprocess==1
if verbosity > 1
disp('************************** Preprocessor ****************************')
end
normpLM=1;
while (k < presteps) && (Psix > eps2) && (normpLM>null)
k = k+1;
% choice of Levenberg-Marquardt parameter, note that we do not use
% the condition estimator for large-scale problems, although this
% may cause numerical problems in some examples
i = false;
mu = 0;
if n<100
i = true;
mu = 1e-16;
if condest(DPhix'*DPhix)>1e25
mu = 1e-6/(k+1);
end
end
if i
pLM = [DPhix; sqrt(mu)*speye(n)]\[-Phix; zeros(n,1)];
else
pLM = -DPhix\Phix;
end
normpLM = norm(pLM);
% compute the projected Levenberg-Marquard step onto box Xk
lbnew = max(min(lb-x,0),-delta);
ubnew = min(max(ub-x,0),delta);
d = max(lbnew,min(pLM,ubnew));
xnew = x+d;
% function evaluations etc.
[Fxnew,DFxnew] = feval(FUN,xnew,varargin{:});
Phixnew = Phi(xnew,Fxnew,lb,ub,lambda1,lambda2,n,Indexset);
Psixnew = 0.5*(Phixnew'*Phixnew);
normPhixnew = norm(Phixnew);
% update of delta
if normPhixnew<=eta*normPhix
delta = max(deltamin,sigma2*delta);
elseif normPhixnew>5*eta*normPhix
delta = max(deltamin,sigma1*delta);
end
% update
x = xnew;
Fx = Fxnew;
DFx = DFxnew;
Phix = Phixnew;
Psix = Psixnew;
normPhix = normPhixnew;
DPhix = DPhi(x,Fx,DFx,lb,ub,lambda1,lambda2,n,Indexset);
DPsix = DPhix'*Phix;
normDPsix = norm(DPsix,inf);
% output at each iteration
t=1;
if verbosity > 1
fprintf('%4.0f %24.5e %24.5e %11.7g\n',k,Psix,normDPsix,t);
end
end
end
% terminate program or redefine current iterate as original initial point
if preprocess==1 && Psix<eps2
if verbosity > 0
fprintf('Psix = %1.4e\nnormDPsix = %1.4e\n',Psix,normDPsix);
disp('Approximate solution found.')
end
EXITFLAG = 1;
FVAL = Fx;
OUTPUT.iterations = k;
OUTPUT.Psix = Psix;
OUTPUT.normDPsix = normDPsix;
JACOB = DFx;
return
elseif preprocess==1 && Psix>=eps2
x=x0;
Phix=Phix0;
Psix=Psix0;
DPhix=DPhix0;
DPsix=DPsix0;
if verbosity > 1
disp('******************** Restart with initial point ********************')
fprintf('%4.0f %24.5e %24.5e\n',k_main,Psix0,normDPsix0);
end
end
%% Main algorithm
if verbosity > 1
disp('************************** Main program ****************************')
end
while (k < kmax) && (Psix > eps2)
% choice of Levenberg-Marquardt parameter, note that we do not use
% the condition estimator for large-scale problems, although this
% may cause numerical problems in some examples
i = false;
if n<100
i = true;
mu = 1e-16;
if condest(DPhix'*DPhix)>1e25
mu = 1e-1/(k+1);
end
end
% compute a Levenberg-Marquard direction
if i
d = [DPhix; sqrt(mu)*speye(n)]\[-Phix; zeros(n,1)];
else
d = -DPhix\Phix;
end
% computation of steplength t using the nonmonotone Armijo-rule
% starting with the 6-th iteration
% computation of steplength t using the monotone Armijo-rule if
% d is a 'good' descent direction or k<=5
t = 1;
xnew = x+d;
Fxnew = feval(FUN,xnew,varargin{:});
Phixnew = Phi(xnew,Fxnew,lb,ub,lambda1,lambda2,n,Indexset);
Psixnew = 0.5*(Phixnew'*Phixnew);
const = sigma*DPsix'*d;
while (Psixnew > MaxPsi + const*t) && (t > tmin)
t = t*beta;
xnew = x+t*d;
Fxnew = feval(FUN,xnew,varargin{:});
Phixnew = Phi(xnew,Fxnew,lb,ub,lambda1,lambda2,n,Indexset);
Psixnew = 0.5*(Phixnew'*Phixnew);
end
% updatings
x = xnew;
Fx = Fxnew;
Phix = Phixnew;
Psix = Psixnew;
[~,DFx] = feval(FUN,x,varargin{:});
DPhix = DPhi(x,Fx,DFx,lb,ub,lambda1,lambda2,n,Indexset);
DPsix = DPhix'*Phix;
normDPsix = norm(DPsix);
k = k+1;
k_main = k_main+1;
if k_main<=5
aux(mod(k_main,m)+1) = Psix;
MaxPsi = Psix;
else
aux(mod(k_main,m)+1) = Psix;
MaxPsi = max(aux);
end
% updatings for the watchdog strategy
if watchdog ==1
if Psix<Psibest
kbest = k;
xbest = x;
Phibest = Phix;
Psibest = Psix;
DPhibest = DPhix;
DPsibest = DPsix;
normDPsibest = normDPsix;
elseif k-kbest>kwatch
x=xbest;
Phix=Phibest;
Psix=Psibest;
DPhix=DPhibest;
DPsix=DPsibest;
normDPsix=normDPsibest;
MaxPsi=Psix;
end
end
if verbosity > 1
% output at each iteration
fprintf('%4.0f %24.5e %24.5e %11.7g\n',k,Psix,normDPsix,t);
end
end
%% Final output
if Psix<=eps2
EXITFLAG = 1;
if verbosity > 0, disp('Approximate solution found.'); end
elseif k>=kmax
EXITFLAG = 0;
if verbosity > 0, disp('Maximum iteration number reached.'); end
elseif normDPsix<=eps1
EXITFLAG = -1; % Provisoire
if verbosity > 0, disp('Approximate stationary point found.'); end
else
EXITFLAG = -1; % Provisoire
if verbosity > 0, disp('No solution found.'); end
end
FVAL = Fx;
OUTPUT.iterations = k;
OUTPUT.Psix = Psix;
OUTPUT.normDPsix = normDPsix;
JACOB = DFx;
%% Subfunctions
function y = Phi(x,Fx,lb,ub,lambda1,lambda2,n,Indexset)
%% PHI
y = zeros(2*n,1);
phi_u = sqrt((ub-x).^2+Fx.^2)-ub+x+Fx;
LZ = false(n,1); % logical zero
I0 = Indexset==0;
y(I0) = -lambda1*Fx(I0);
y([LZ; I0]) = -lambda2*Fx(I0);
I1 = Indexset==1;
y(I1) = lambda1*(-x(I1)+lb(I1)-Fx(I1)+sqrt((x(I1)-lb(I1)).^2+Fx(I1).^2));
y([LZ; I1]) = lambda2*max(0,x(I1)-lb(I1)).*max(0,Fx(I1));
I2 = Indexset==2;
y(I2) = -lambda1*phi_u(I2);
y([LZ; I2]) = lambda2*max(0,ub(I2)-x(I2)).*max(0,-Fx(I2));
I3 = Indexset==3;
y(I3) = lambda1*(sqrt((x(I3)-lb(I3)).^2+phi_u(I3).^2)-x(I3)+lb(I3)-phi_u(I3));
y([LZ; I3]) = lambda2*(max(0,x(I3)-lb(I3)).*max(0,Fx(I3))+max(0,ub(I3)-x(I3)).*max(0,-Fx(I3)));
function H = DPhi(x,Fx,DFx,lb,ub,lambda1,lambda2,n,Indexset)
%% DPHI evaluates an element of the C-subdifferential of operator Phi
null = 1e-8;
beta_l = zeros(n,1);
beta_u = zeros(n,1);
alpha_l = zeros(n,1);
alpha_u = zeros(n,1);
z = zeros(n,1);
H2 = sparse(n,n);
I = abs(x-lb)<=null & abs(Fx)<=null;
beta_l(I) = 1;
z(I) = 1;
I = abs(ub-x)<=null & abs(Fx)<=null;
beta_u(I) = 1;
z(I) = 1;
I = x-lb>=-null & Fx>=-null;
alpha_l(I) = 1;
I = ub-x>=-null & Fx<=null;
alpha_u(I) = 1;
Da = zeros(n,1);
Db = zeros(n,1);
I = 1:n;
I0 = Indexset==0;
Da(I0) = 0;
Db(I0) = -1;
H2(I0,:) = -DFx(I0,:);
I1 = Indexset==1;
denom1 = zeros(n,1);
denom2 = zeros(n,1);
if any(I1)
denom1(I1) = max(null,sqrt((x(I1)-lb(I1)).^2+Fx(I1).^2));
denom2(I1) = max(null,sqrt(z(I1).^2+(DFx(I1,:)*z).^2));
end
I1b = Indexset==1 & beta_l==0;
Da(I1b) = (x(I1b)-lb(I1b))./denom1(I1b)-1;
Db(I1b) = Fx(I1b)./denom1(I1b)-1;
I1b = Indexset==1 & beta_l~=0;
if any(I1b)
Da(I1b) = z(I1b)./denom2(I1b)-1;
Db(I1b) = (DFx(I1b,:)*z)./denom2(I1b)-1;
end
I1a = I(Indexset==1 & alpha_l==1);
if any(I1a)
H2(I1a,:) = bsxfun(@times,x(I1a)-lb(I1a),DFx(I1a,:))+...
sparse(1:length(I1a),I1a,Fx(I1a),length(I1a),n,length(I1a));
end
I2 = Indexset==2;
denom1 = zeros(n,1);
denom2 = zeros(n,1);
if any(I2)
denom1(I2) = max(null,sqrt((ub(I2)-x(I2)).^2+Fx(I2).^2));
denom2(I2) = max(null,sqrt(z(I2).^2+(DFx(I2,:)*z).^2));
end
I2b = Indexset==2 & beta_u==0;
Da(I2b) = (ub(I2b)-x(I2b))./denom1(I2b)-1;
Db(I2b) = -Fx(I2b)./denom1(I2b)-1;
I2b = Indexset==2 & beta_u~=0;
if any(I2b)
Da(I2b) = -z(I2b)./denom2(I2b)-1;
Db(I2b) = -(DFx(I2b,:)*z)./denom2(I2b)-1;
end
I2a = I(Indexset==2 & alpha_u==1);
if any(I2a)
H2(I2a,:) = bsxfun(@times,x(I2a)-ub(I2a),DFx(I2a,:))+...
sparse(1:length(I2a),I2a,Fx(I2a),length(I2a),n,length(I2a));
end
I3 = Indexset==3;
phi = zeros(n,1);
ai = zeros(n,1);
bi = zeros(n,1);
ci = zeros(n,1);
di = zeros(n,1);
denom1 = zeros(n,1);
denom2 = zeros(n,1);
denom3 = zeros(n,1);
denom4 = zeros(n,1);
if any(I3)
phi(I3) = -ub(I3)+x(I3)+Fx(I3)+sqrt((ub(I3)-x(I3)).^2+Fx(I3).^2);
denom1(I3) = max(null,sqrt((x(I3)-lb(I3)).^2+phi(I3).^2));
denom2(I3) = max(null,sqrt(z(I3).^2+(DFx(I3,:)*z).^2));
denom3(I3) = max(null,sqrt((ub(I3)-x(I3)).^2+Fx(I3).^2));
denom4(I3) = max(null,sqrt(z(I3).^2));
end
I3bu = Indexset==3 & beta_u==0;
ci(I3bu) = (x(I3bu)-ub(I3bu))./denom3(I3bu)+1;
di(I3bu) = Fx(I3bu)./denom3(I3bu)+1;
I3bu = Indexset==3 & beta_u~=0;
if any(I3bu)
ci(I3bu) = 1+z(I3bu)./denom2(I3bu);
di(I3bu) = 1+(DFx(I3bu,:)*z)./denom2(I3bu);
end
I3bl = Indexset==3 & beta_l==0;
ai(I3bl) = (x(I3bl)-lb(I3bl))./denom1(I3bl)-1;
bi(I3bl) = phi(I3bl)./denom1(I3bl)-1;
I3bl = Indexset==3 & beta_l~=0;
if any(I3bl)
ai(I3bl) = z(I3bl)./denom4(I3bl)-1;
bi(I3bl) = (ci(I3bl).*z(I3bl)+(di(I3bl,ones(1,n)).*DFx(I3bl,:))*z)./denom4(I3bl)-1;
end
Da(I3) = ai(I3)+bi(I3).*ci(I3);
Db(I3) = bi(I3).*di(I3);
I3a = I(Indexset==3 & alpha_l==1 & alpha_u==1);
if any(I3a)
H2(I3a,:) = bsxfun(@times,-lb(I3a)-ub(I3a)+2*x(I3a),DFx(I3a,:))+...
2*sparse(1:length(I3a),I3a,Fx(I3a),length(I3a),n,length(I3a));
end
I3a = I(Indexset==3 & alpha_l==1 & alpha_u~=1);
if any(I3a)
H2(I3a,:) = bsxfun(@times,x(I3a)-lb(I3a),DFx(I3a,:))+...
sparse(1:length(I3a),I3a,Fx(I3a),length(I3a),n,length(I3a));
end
I3a = I(Indexset==3 & alpha_l~=1 & alpha_u==1);
if any(I3a)
H2(I3a,:) = bsxfun(@times,x(I3a)-ub(I3a),DFx(I3a,:))+...
sparse(1:length(I3a),I3a,Fx(I3a),length(I3a),n,length(I3a));
end
H1 = bsxfun(@times,Db,DFx);
H1 = spdiags(diag(H1)+Da,0,H1);
H = [lambda1*H1; lambda2*H2];

71
matlab/perfect_foresight_problem.m Normal file
View File

@ -0,0 +1,71 @@
function [residuals,JJacobian] = perfect_foresight_problem(y, dynamic_function, Y0, YT, ...
exo_simul, params, steady_state, ...
T, ny, i_cols, ...
i_cols_J1, i_cols_1, i_cols_T, ...
i_cols_j,nnzJ)
% function perfect_foresight_problem(x, model_dynamic, Y0, YT,exo_simul,
% params, steady_state, periods, ny, i_cols,i_cols_J1, i_cols_1,
% i_cols_T, i_cols_j, nnzA)
% computes the residuals and th Jacobian matrix
% for a perfect foresight problem over T periods.
%
% INPUTS
% ...
% OUTPUTS
% ...
% ALGORITHM
% ...
%
% SPECIAL REQUIREMENTS
% None.
% Copyright (C) 1996-2014 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
YY = [Y0; y; YT];
residuals = zeros(T*ny,1);
if nargout == 2
JJacobian = sparse([],[],[],T*ny,T*ny,T*nnzJ);
end
i_rows = 1:ny;
i_cols_J = i_cols;
for it = 2:(T+1)
if nargout == 1
residuals(i_rows) = dynamic_function(YY(i_cols),exo_simul, params, ...
steady_state,it);
elseif nargout == 2
[residuals(i_rows),jacobian] = dynamic_function(YY(i_cols),exo_simul, params, ...
steady_state,it);
if it == 2
JJacobian(i_rows,i_cols_J1) = jacobian(:,i_cols_1);
elseif it == T + 1
JJacobian(i_rows,i_cols_J(i_cols_T)) = jacobian(:,i_cols_T);
else
JJacobian(i_rows,i_cols_J) = jacobian(:,i_cols_j);
i_cols_J = i_cols_J + ny;
end
end
i_rows = i_rows + ny;
i_cols = i_cols + ny;
end

37
matlab/perfect_foresight_solver.m
View File

@ -31,12 +31,12 @@ function perfect_foresight_solver()
global M_ options_ oo_
if options_.stack_solve_algo < 0 || options_.stack_solve_algo > 6
error('PERFECT_FORESIGHT_SOLVER: stack_solve_algo must be between 0 and 6')
if options_.stack_solve_algo < 0 || options_.stack_solve_algo > 7
error('PERFECT_FORESIGHT_SOLVER: stack_solve_algo must be between 0 and 7')
end
if ~options_.block && ~options_.bytecode && options_.stack_solve_algo ~= 0 ...
&& options_.stack_solve_algo ~= 6
&& options_.stack_solve_algo ~= 6 && options_.stack_solve_algo ~= 7
error('PERFECT_FORESIGHT_SOLVER: you must use stack_solve_algo=0 or stack_solve_algo=6 when not using block nor bytecode option')
end
@ -185,8 +185,37 @@ else
else % General case
if options_.stack_solve_algo == 0
sim1;
else % stack_solve_algo = 6
elseif options_.stack_solve_algo == 6
sim1_lbj;
elseif options_.stack_solve_algo == 7
periods = options_.periods;
if ~isfield(options_.lmmcp,'lb')
[lb,ub,pfm.eq_index] = get_complementarity_conditions(M_);
options_.lmmcp.lb = repmat(lb,periods,1);
options_.lmmcp.ub = repmat(ub,periods,1);
end
y = oo_.endo_simul;
y0 = y(:,1);
yT = y(:,periods+2);
z = y(:,2:periods+1);
illi = M_.lead_lag_incidence';
[i_cols,~,i_cols_j] = find(illi(:));
illi = illi(:,2:3);
[i_cols_J1,~,i_cols_1] = find(illi(:));
i_cols_T = nonzeros(M_.lead_lag_incidence(1:2,:)');
[y,info] = dynare_solve(@perfect_foresight_problem,z(:),1, ...
str2func([M_.fname '_dynamic']),y0,yT, ...
oo_.exo_simul,M_.params,oo_.steady_state, ...
options_.periods,M_.endo_nbr,i_cols, ...
i_cols_J1, i_cols_1, i_cols_T, i_cols_j, ...
M_.NNZDerivatives(1));
oo_.endo_simul = [y0 reshape(y,M_.endo_nbr,periods) yT];
if info == 1
oo_.deterministic_simulation.status = 0;
else
oo_.deterministic_simulation.status = 1;
end;
end
end
end

13
matlab/reports/@graph/graph.m
View File

@ -72,6 +72,9 @@ o.xTickLabelAnchor = 'east';
o.width = 6;
o.height = 4.5;
o.miscTikzPictureOptions = '';
o.miscTikzAxisOptions = '';
if nargin == 1
assert(isa(varargin{1}, 'graph'),['@graph.graph: with one arg you ' ...
'must pass a graph object']);
@ -87,7 +90,7 @@ elseif nargin > 1
% overwrite default values
for pair = reshape(varargin, 2, [])
ind = strmatch(lower(pair{1}), lower(optNames), 'exact');
ind = find(strcmpi(optNames, pair{1}));
assert(isempty(ind) || length(ind) == 1);
if ~isempty(ind)
o.(optNames{ind}) = pair{2};
@ -105,6 +108,8 @@ assert(iscellstr(o.title), '@graph.graph: title must be a cell array of string(s
assert(ischar(o.titleFormat), '@graph.graph: titleFormat file must be a string');
assert(ischar(o.xlabel), '@graph.graph: xlabel file must be a string');
assert(ischar(o.ylabel), '@graph.graph: ylabel file must be a string');
assert(ischar(o.miscTikzPictureOptions), '@graph.graph: miscTikzPictureOptions file must be a string');
assert(ischar(o.miscTikzAxisOptions), '@graph.graph: miscTikzAxisOptions file must be a string');
assert(ischar(o.graphName), '@graph.graph: graphName must be a string');
assert(ischar(o.graphDirName), '@graph.graph: graphDirName must be a string');
assert(islogical(o.showGrid), '@graph.graph: showGrid must be either true or false');
@ -141,17 +146,17 @@ valid_legend_orientations = {'vertical', 'horizontal'};
assert(any(strcmp(o.legendOrientation, valid_legend_orientations)), ...
['@graph.graph: legendOrientation must be one of ' strjoin(valid_legend_orientations, ' ')]);
assert(isempty(o.shade) || (isa(o.shade, 'dates') && o.shade.ndat >= 2), ...
assert(isempty(o.shade) || (isdates(o.shade) && o.shade.ndat >= 2), ...
['@graph.graph: shade is specified as a dates range, e.g. ' ...
'''dates(''1999q1''):dates(''1999q3'')''.']);
assert(isempty(o.xrange) || (isa(o.xrange, 'dates') && o.xrange.ndat >= 2), ...
assert(isempty(o.xrange) || (isdates(o.xrange) && o.xrange.ndat >= 2), ...
['@graph.graph: xrange is specified as a dates range, e.g. ' ...
'''dates(''1999q1''):dates(''1999q3'')''.']);
assert(isempty(o.yrange) || (isfloat(o.yrange) && length(o.yrange) == 2 && ...
o.yrange(1) < o.yrange(2)), ...
['@graph.graph: yrange is specified an array with two float entries, ' ...
'the lower bound and upper bound.']);
assert(isempty(o.data) || isa(o.data, 'dseries'), ['@graph.graph: data must ' ...
assert(isempty(o.data) || isdseries(o.data), ['@graph.graph: data must ' ...
'be a dseries']);
assert(isempty(o.seriesToUse) || iscellstr(o.seriesToUse), ['@graph.graph: ' ...
'seriesToUse must be a cell array of string(s)']);

22
matlab/reports/@graph/writeGraphFile.m
View File

@ -49,7 +49,11 @@ if fid == -1
error(['@graph.writeGraphFile: ' msg]);
end
fprintf(fid, '\\begin{tikzpicture}[baseline]');
fprintf(fid, '\\begin{tikzpicture}[baseline');
if ~isempty(o.miscTikzPictureOptions)
fprintf(fid, ',%s', o.miscTikzPictureOptions);
end
fprintf(fid, ']');
if isempty(o.xrange)
dd = getMaxRange(o.series);
@ -141,6 +145,7 @@ if isempty(o.yrange)
else
fprintf(fid, 'ymin=%f,\nymax=%f,\n',o.yrange(1),o.yrange(2));
end
fprintf(fid, 'xmin = 1,\nxmax = %d,\n', length(dd));
if o.showLegend
fprintf(fid, 'legend style={');
@ -165,6 +170,10 @@ end
if ~isempty(o.ylabel)
fprintf(fid, 'ylabel=%s,\n', o.ylabel);
end
if ~isempty(o.miscTikzAxisOptions)
fprintf(fid, '%s', o.miscTikzAxisOptions);
end
fprintf(fid, ']\n');
if o.showZeroline
@ -174,7 +183,10 @@ end
for i=1:ne
o.series{i}.writeSeriesForGraph(fid, dd);
if o.showLegend
fprintf(fid, '\\addlegendentry{%s}\n', o.series{i}.getTexName());
le = o.series{i}.getNameForLegend();
if ~isempty(le)
fprintf(fid, '\\addlegendentry{%s}\n', le);
end
end
end
@ -186,18 +198,18 @@ if ~isempty(o.shade)
date2string(o.shade(1)) ' or ' date2string(o.shade(end)) ' is not in the date ' ...
'range of data selected.']);
if x1 == 1
fprintf(fid,['\\begin{pgfonlayer}{background}\n\\fill[%s!%f]\n(axis ' ...
fprintf(fid,['\\begin{pgfonlayer}{background0}\n\\fill[%s!%f]\n(axis ' ...
'cs:\\pgfkeysvalueof{/pgfplots/xmin},\\pgfkeysvalueof{/pgfplots/ymin})\nrectangle (axis ' ...
'cs:%f,\\pgfkeysvalueof{/pgfplots/ymax});\n\\end{pgfonlayer}\n'], ...
o.shadeColor, o.shadeOpacity, x2);
elseif x2 == dd.ndat
fprintf(fid,['\\begin{pgfonlayer}{background}\n\\fill[%s!%f]\n(axis ' ...
fprintf(fid,['\\begin{pgfonlayer}{background0}\n\\fill[%s!%f]\n(axis ' ...
'cs:%f,\\pgfkeysvalueof{/pgfplots/ymin})\nrectangle (axis ' ...
'cs:\\pgfkeysvalueof{/pgfplots/xmax},\\pgfkeysvalueof{/' ...
'pgfplots/ymax});\n\\end{pgfonlayer}\n'], ...
o.shadeColor, o.shadeOpacity, x1);
else
fprintf(fid,['\\begin{pgfonlayer}{background}\n\\fill[%s!%f]\n(axis ' ...
fprintf(fid,['\\begin{pgfonlayer}{background0}\n\\fill[%s!%f]\n(axis ' ...
'cs:%f,\\pgfkeysvalueof{/pgfplots/ymin})\nrectangle (axis ' ...
'cs:%f,\\pgfkeysvalueof{/pgfplots/ymax});\n\\end{pgfonlayer}\n'], ...
o.shadeColor, o.shadeOpacity, x1, x2);

2
matlab/reports/@page/page.m
View File

@ -53,7 +53,7 @@ elseif nargin > 1
% overwrite default values
for pair = reshape(varargin, 2, [])
ind = strmatch(lower(pair{1}), lower(optNames), 'exact');
ind = find(strcmpi(optNames, pair{1}));
assert(isempty(ind) || length(ind) == 1);
if ~isempty(ind)
o.(optNames{ind}) = pair{2};

32
matlab/reports/@paragraph/display.m Normal file
View File

@ -0,0 +1,32 @@
function display(o)
%function display(o)
% Display a Paragraph object
%
% INPUTS
% o [paragraph] paragraph object
%
% OUTPUTS
% none
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2014 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
display_reporting_object(o);
end

63
matlab/reports/@paragraph/paragraph.m Normal file
View File

@ -0,0 +1,63 @@
function o = paragraph(varargin)
%function o = paragraph(varargin)
% Instantiates a paragraph object
% Copyright (C) 2014 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
o = struct;
o.balancedCols = false;
o.cols = 1;
o.heading = '';
o.indent = true;
o.text = '';
if nargin == 1
assert(isa(varargin{1}, 'paragraph'),['With one arg to Paragraph constructor, ' ...
'you must pass a paragraph object']);
o = varargin{1};
return;
elseif nargin > 1
if round(nargin/2) ~= nargin/2
error(['@paragraph.paragraph: options must be supplied in name/value ' ...
'pairs.']);
end
optNames = fieldnames(o);
% overwrite default values
for pair = reshape(varargin, 2, [])
ind = find(strcmpi(optNames, pair{1}));
assert(isempty(ind) || length(ind) == 1);
if ~isempty(ind)
o.(optNames{ind}) = pair{2};
else
error('@paragraph.paragraph: %s is not a recognized option.', pair{1});
end
end
end
assert(islogical(o.indent), '@paragraph.paragraph: indent must be either true or false');
assert(islogical(o.balancedCols), '@paragraph.paragraph: balancedCols must be either true or false');
assert(isint(o.cols), '@paragraph.paragraph: cols must be an integer');
assert(ischar(o.text), '@paragraph.paragraph: text must be a string');
assert(ischar(o.heading), '@paragraph.paragraph: heading must be a string');
% Create paragraph object
o = class(o, 'paragraph');
end

50
matlab/reports/@paragraph/subsasgn.m Normal file
View File

@ -0,0 +1,50 @@
function B = subsasgn(A, S, V)
% function B = subsasgn(A, S, V)
% Copyright (C) 2014 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
B = A;
if length(S) > 1
for i=1:(length(S)-1)
B = subsref(B, S(i));
end
B = subsasgn(B, S(end), V);
B = subsasgn(A, S(1:(end-1)), B);
return
end
switch S.type
case '()'
index = S.subs{:};
assert(isnumeric(index));
B.elements{index} = V;
case '{}'
index = S.subs{:};
assert(isnumeric(index));
B{index} = V;
case '.'
switch S.subs
case fieldnames(A)
B.(S.subs) = V;
otherwise
error(['@paragraph.subsasgn: field ' S.subs 'does not exist']);
end
otherwise
error('@paragraph.subsasgn: syntax error');
end
end

48
matlab/reports/@paragraph/subsref.m Normal file
View File

@ -0,0 +1,48 @@
function A = subsref(A, S)
%function A = subsref(A, S)
% Copyright (C) 2014 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
switch S(1).type
case '.'
switch S(1).subs
case fieldnames(A)
A = A.(S(1).subs);
case methods(A)
if areParensNext(S)
A = feval(S(1).subs, A, S(2).subs{:});
S = shiftS(S,1);
else
A = feval(S(1).subs, A);
end
otherwise
error(['@paragraph.subsref: unknown field or method: ' S(1).subs]);
end
case '()'
A = A.elements{S(1).subs{:}};
case '{}'
error(['@paragraph.subsref: ' S(1).type ' indexing not supported.']);
otherwise
error('@paragraph.subsref: impossible case')
end
S = shiftS(S,1);
if length(S) >= 1
A = subsref(A, S);
end
end

70
matlab/reports/@paragraph/write.m Normal file
View File

@ -0,0 +1,70 @@
function o = write(o, fid)
%function o = write(o, fid)
% Write Paragraph object
%
% INPUTS
% o [paragraph] paragraph object
% fid [integer] file id
%
% OUTPUTS
% o [paragraph] paragraph object
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2014 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
assert(fid ~= -1);
fprintf(fid, '%% Paragraph Object\n\\multicolumn{1}{p{\\textwidth}}{%%\n');
if o.cols ~= 1
bc = '';
if o.balancedCols
bc = '*';
end
fprintf(fid, '\\begin{multicols%s}{%d}%%\n', bc, o.cols);
end
if ~isempty(o.heading)
if o.cols ~= 1
fprintf(fid, '[%s\n]\n', o.heading);
else
fprintf(fid, '%s\\newline \\newline\n', o.heading);
end
end
if o.indent
fprintf(fid, '\\hspace{4ex}');
end
fprintf(fid, '%s', o.text);
if o.cols ~= 1
fprintf(fid, '\\end{multicols%s}\n', bc);
end
fprintf(fid, '}\n%% End Paragraph Object\n\n');
end
%\multicolumn{1}{p{\textwidth}}
%{\begin{multicols}{2}
%Hello, here is some text without a meaning. This text should show what
%a printed text will look like at this place.
%\columnbreak
%If you read this text, you will get no information. Really? Is there
%no information? Is there...
%\end{multicols}}\\

34
matlab/reports/@report/addParagraph.m Normal file
View File

@ -0,0 +1,34 @@
function o = addParagraph(o, varargin)
%function o = addParagraph(o, varargin)
% Add a paragraph to the current section of the current page in the report
%
% INPUTS
% o [report] report object
% varargin arguments to @section/addGraph.m
%
% OUTPUTS
% o [report] updated report object
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2013-2014 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
o.pages{end}.sections{end} = ...
o.pages{end}.sections{end}.addParagraph(varargin{:});
end

2
matlab/reports/@report/compile.m
View File

@ -79,7 +79,7 @@ if status ~= 0
error(['@report.compile: There was an error in compiling ' rfn '.pdf.' ...
' ' compiler ' returned the error code: ' num2str(status)]);
end
fprintf(1, '\n\nDone.\n')
fprintf(1, '\n\nDone.\n');
disp('Your compiled report is located here:');
disp([' ' pwd filesep rfn '.pdf']);

2
matlab/reports/@report/report.m
View File

@ -57,7 +57,7 @@ elseif nargin > 1
% overwrite default values
for pair = reshape(varargin, 2, [])
ind = strmatch(lower(pair{1}), lower(optNames), 'exact');
ind = find(strcmpi(optNames, pair{1}));
assert(isempty(ind) || length(ind) == 1);
if ~isempty(ind)
o.(optNames{ind}) = pair{2};

8
matlab/reports/@report/write.m
View File

@ -41,7 +41,7 @@ if strcmpi(o.orientation, 'landscape')
fprintf(fid, ',landscape');
end
fprintf(fid, ']{geometry}\n');
fprintf(fid, '\\usepackage{pdflscape, booktabs, pgfplots, colortbl, adjustbox}\n');
fprintf(fid, '\\usepackage{pdflscape, booktabs, pgfplots, colortbl, adjustbox, multicol}\n');
fprintf(fid, '\\pgfplotsset{compat=1.5.1}');
fprintf(fid, ['\\makeatletter\n' ...
'\\def\\blfootnote{\\gdef\\@thefnmark{}\\@footnotetext}\n' ...
@ -71,9 +71,9 @@ fprintf(fid, '\\renewcommand{\\bottomfraction}{0.8}\n');
fprintf(fid, '\\setlength{\\parindent}{0in}\n');
fprintf(fid, '\\newlength\\sectionheight\n');
fprintf(fid, '\\begin{document}\n');
fprintf(fid, '\\pgfdeclarelayer{background}\n');
fprintf(fid, '\\pgfdeclarelayer{foreground}\n');
fprintf(fid, '\\pgfsetlayers{background,main,foreground}\n');
fprintf(fid, '\\pgfdeclarelayer{background0}\n');
fprintf(fid, '\\pgfdeclarelayer{background1}\n');
fprintf(fid, '\\pgfsetlayers{background0,background1,main}\n');
fprintf(fid, '\\centering\n');
nps = length(o.pages);

31
matlab/reports/@report_series/getNameForLegend.m Normal file
View File

@ -0,0 +1,31 @@
function s = getNameForLegend(o)
%function s = getNameForLegend(o)
% Copyright (C) 2014 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if isempty(o.data) || ~o.graphShowInLegend
% for the case when there is no data in the series
% e.g. graphVline was passed
% or when the user does not want this series shown in
% the legend
s = '';
else
assert(size(o.data,2) == 1);
s = o.data.tex{:};
end
end

8
matlab/reports/@report_series/getRange.m
View File

@ -18,6 +18,10 @@ function dd = getRange(o)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
assert(~isempty(o.data) && size(o.data, 2) == 1);
dd = o.data.dates;
if isempty(o.data)
dd = dates();
else
assert(size(o.data, 2) == 1);
dd = o.data.dates;
end
end

12
matlab/reports/@report_series/getTexName.m
View File

@ -1,7 +1,7 @@
function s = getTexName(o)
%function s = getTexName(o)
% Copyright (C) 2013 Dynare Team
% Copyright (C) 2013-2014 Dynare Team
%
% This file is part of Dynare.
%
@ -18,6 +18,12 @@ function s = getTexName(o)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
assert(~isempty(o.data) && size(o.data, 2) == 1);
s = o.data.tex{:};
if isempty(o.data)
% for the case when there is no data in the series
% e.g. graphVline was passed
s = '';
else
assert(size(o.data,2) == 1);
s = o.data.tex{:};
end
end

2
matlab/reports/@report_series/printSeries.m
View File

@ -36,7 +36,7 @@ dataString = sprintf('%%.%df', precision);
precision = 10^precision;
data = dser(dates);
data = data.data;
data = setDataToZeroFromZeroTol(o, data);
for i=1:size(data,1)
fprintf(fid, '&');
if o.tableShowMarkers

15
matlab/reports/@report_series/report_series.m
View File

@ -35,15 +35,24 @@ o = struct;
o.data = '';
o.graphLegendName = '';
o.graphLineColor = 'black';
o.graphLineStyle = 'solid';
o.graphLineWidth = 0.5;
o.graphShowInLegend = true;
o.graphMarker = '';
o.graphMarkerEdgeColor = '';
o.graphMarkerFaceColor = '';
o.graphMarkerSize = 1;
o.graphMiscTikzAddPlotOptions = '';
o.graphHline = {};
o.graphVline = dates();
o.tableShowMarkers = false;
o.tableNegColor = 'red';
o.tablePosColor = 'blue';
@ -74,7 +83,7 @@ elseif nargin > 1
% overwrite default values
for pair = reshape(varargin, 2, [])
ind = strmatch(lower(pair{1}), lower(optNames), 'exact');
ind = find(strcmpi(optNames, pair{1}));
assert(isempty(ind) || length(ind) == 1);
if ~isempty(ind)
o.(optNames{ind}) = pair{2};
@ -84,6 +93,10 @@ elseif nargin > 1
end
end
if ~isempty(o.graphLegendName)
o.data = o.data.tex_rename(o.graphLegendName);
end
% Create report_series object
o = class(o, 'report_series');
end

28
matlab/reports/@report_series/setDataToZeroFromZeroTol.m Normal file
View File

@ -0,0 +1,28 @@
function d = setDataToZeroFromZeroTol(o, ds)
%function d = setDataToZeroFromZeroTol(o, ds)
% Copyright (C) 2014 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
d = ds.data;
stz = bsxfun(@and, ...
bsxfun(@lt, d, o.zeroTol), ...
bsxfun(@gt, d, -o.zeroTol));
if any(stz)
d(stz) = 0;
end
end

70
matlab/reports/@report_series/writeSeriesForGraph.m
View File

@ -30,8 +30,14 @@ function o = writeSeriesForGraph(o, fid, xrange)
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
%% Validate options provided by user
assert(~isempty(o.data) && isa(o.data, 'dseries'), ['@report_series.writeSeriesForGraph: must ' ...
'provide data as a dseries']);
if isempty(o.graphVline) && isempty(o.graphHline)
assert(~isempty(o.data) && isdseries(o.data), ['@report_series.writeSeriesForGraph: must ' ...
'provide data as a dseries']);
end
assert(ischar(o.graphMiscTikzAddPlotOptions), ['@report_series.writeSeriesForGraph: ' ...
'graphMiscTikzAddPlotOptions file must be a string']);
assert(islogical(o.graphShowInLegend), '@graph.graph: graphShowInLegend must be either true or false');
% Line
valid_graphLineColor = {'red', 'green', 'blue', 'cyan ', 'magenta', 'yellow', ...
@ -61,12 +67,37 @@ assert(isfloat(o.graphMarkerSize) && o.graphMarkerSize > 0, ...
assert(~(strcmp(o.graphLineStyle, 'none') && isempty(o.graphMarker)), ['@report_series.writeSeriesForGraph: ' ...
'you must provide at least one of graphLineStyle and graphMarker']);
% Validate xrange
% Validate graphVline
assert(isempty(o.graphVline) || (isdates(o.graphVline) && o.graphVline.ndat == 1), ...
'@report_series.writeSeriesForGraph: graphVline must be a dates of size one');
assert(isempty(o.graphHline) || isnumeric(o.graphHline), ...
'@report_series.writeSeriesForGraph: graphHline must a single numeric value');
% Zero tolerance
assert(isfloat(o.zeroTol), '@report_series.write: zeroTol must be a float');
%% graphVline && graphHline
if ~isempty(o.graphVline)
fprintf(fid, '%%Vertical Line\n\\begin{pgfonlayer}{background1}\n\\draw');
writeLineOptions(o, fid);
stringsdd = strings(xrange);
x = find(strcmpi(date2string(o.graphVline), stringsdd));
fprintf(fid, ['(axis cs:%d,\\pgfkeysvalueof{/pgfplots/ymin}) -- (axis ' ...
'cs:%d,\\pgfkeysvalueof{/pgfplots/ymax});\n\\end{pgfonlayer}\n'], ...
x, x);
end
if ~isempty(o.graphHline)
fprintf(fid, '%%Horizontal Line\n\\begin{pgfonlayer}{background1}\n\\addplot');
writeLineOptions(o, fid);
fprintf(fid, ['coordinates {(\\pgfkeysvalueof{/pgfplots/xmin},%f)' ...
'(\\pgfkeysvalueof{/pgfplots/xmax},%f)};\n\\end{pgfonlayer}\n'], ...
o.graphHline, o.graphHline);
end
if ~isempty(o.graphVline) || ~isempty(o.graphHline)
% return since the code below assumes that o.data exists
return
end
%%
if isempty(xrange) || all(xrange == o.data.dates)
ds = o.data;
@ -74,18 +105,20 @@ else
ds = o.data(xrange);
end
% if graphing data that is within zeroTol, set to zero, create report_series and
% get line:
thedata = ds.data;
stz = bsxfun(@and, ...
bsxfun(@lt, thedata, o.zeroTol), ...
bsxfun(@gt, thedata, -o.zeroTol));
if any(stz)
thedata(stz) = 0;
thedata = setDataToZeroFromZeroTol(o, ds);
fprintf(fid, '%%series %s\n\\addplot', o.data.name{:});
writeLineOptions(o, fid);
fprintf(fid,'\ntable[row sep=crcr]{\nx y\\\\\n');
for i=1:ds.dates.ndat
fprintf(fid, '%d %f\\\\\n', i, thedata(i));
end
fprintf(fid,'};\n');
end
fprintf(fid, '%%series %s\n\\addplot[color=%s,%s,line width=%fpt,line join=round',...
o.data.name{:}, o.graphLineColor, o.graphLineStyle, o.graphLineWidth);
function writeLineOptions(o, fid)
fprintf(fid, '[color=%s,%s,line width=%fpt,line join=round',...
o.graphLineColor, o.graphLineStyle, o.graphLineWidth);
if ~isempty(o.graphMarker)
if isempty(o.graphMarkerEdgeColor)
o.graphMarkerEdgeColor = o.graphLineColor;
@ -94,11 +127,10 @@ if ~isempty(o.graphMarker)
o.graphMarkerFaceColor = o.graphLineColor;
end
fprintf(fid, ',mark=%s,mark size=%f,every mark/.append style={draw=%s,fill=%s}',...
o.graphMarker,o.graphMarkerSize,o.graphMarkerEdgeColor,o.graphMarkerFaceColor);
o.graphMarker,o.graphMarkerSize,o.graphMarkerEdgeColor,o.graphMarkerFaceColor);
end
fprintf(fid,']\ntable[row sep=crcr]{\nx y\\\\\n');
for i=1:ds.dates.ndat
fprintf(fid, '%d %f\\\\\n', i, thedata(i));
if ~isempty(o.graphMiscTikzAddPlotOptions)
fprintf(fid, ',%s', o.graphMiscTikzAddPlotOptions);
end
fprintf(fid,'};\n');
fprintf(fid,']');
end

6
matlab/reports/@report_series/writeSeriesForTable.m
View File

@ -35,18 +35,18 @@ function o = writeSeriesForTable(o, fid, dates, precision)
%% Validate options passed to function
assert(fid ~= -1);
for i=1:length(dates)
assert(isa(dates{i}, 'dates'));
assert(isdates(dates{i}));
end
assert(isint(precision));
%% Validate options provided by user
assert(ischar(o.tableSubSectionHeader), '@report_series.writeSeriesForTable: tableSubSectionHeader must be a string');
if isempty(o.tableSubSectionHeader)
assert(~isempty(o.data) && isa(o.data, 'dseries'), ...
assert(~isempty(o.data) && isdseries(o.data), ...
'@report_series.writeSeriesForTable: must provide data as a dseries');
if ~isempty(o.tableDataRhs)
assert(~isempty(o.tableDataRhs) && isa(o.tableDataRhs, 'dseries'), ...
assert(~isempty(o.tableDataRhs) && isdseries(o.tableDataRhs), ...
'@report_series.writeSeriesForTable: must provide tableDataRhs as a dseries');
assert(iscell(dates) && length(dates) == 2, ...
'@report_series.writeSeriesForTable: must provide second range with tableDataRhs');

6
matlab/reports/@report_table/report_table.m
View File

@ -65,7 +65,7 @@ elseif nargin > 1
% overwrite default values
for pair = reshape(varargin, 2, [])
ind = strmatch(lower(pair{1}), lower(optNames), 'exact');
ind = find(strcmpi(optNames, pair{1}));
assert(isempty(ind) || length(ind) == 1);
if ~isempty(ind)
o.(optNames{ind}) = pair{2};
@ -78,7 +78,7 @@ if ~iscell(o.range)
o.range = {o.range};
end
if isa(o.vlineAfter, 'dates')
if isdates(o.vlineAfter)
o.vlineAfter = {o.vlineAfter};
end
@ -98,7 +98,7 @@ assert(isint(o.precision), '@report_table.report_table: precision must be an int
assert(isempty(o.range) || length(o.range) <=2 && allCellsAreDatesRange(o.range), ...
['@report_table.report_table: range is specified as a dates range, e.g. ' ...
'''dates(''1999q1''):dates(''1999q3'')''.']);
assert(isempty(o.data) || isa(o.data, 'dseries'), ...
assert(isempty(o.data) || isdseries(o.data), ...
'@report_table.report_table: data must be a dseries');
assert(isempty(o.seriesToUse) || iscellstr(o.seriesToUse), ...
'@report_table.report_table: seriesToUse must be a cell array of string(s)');

4
matlab/reports/@section/addGraph.m
View File

@ -30,5 +30,9 @@ function o = addGraph(o, varargin)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
for i=1:length(o.elements)
assert(~isa(o.elements{i}, 'paragraph'), ...
'@addGraph: A Section that contains a Paragraph cannot contain a Graph');
end
o.elements{end+1} = graph(varargin{:});
end

42
matlab/reports/@section/addParagraph.m Normal file
View File

@ -0,0 +1,42 @@
function o = addParagraph(o, varargin)
%function o = addParagraph(o, varargin)
% Add a paragraph to the Cell Array of elements in this section
%
% INPUTS
% 1 args => add empty paragraph
% 2 args => add given paragraph
% 3 args => add paragraph at index
%
% OUTPUTS
% updated page object
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2014 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
assert(o.cols == 1, ...
['@addParagraph: you can only add a paragraph to a Section that ' ...
'contains one column']);
for i=1:length(o.elements)
assert(isa(o.elements{i}, 'paragraph'), ...
['@addParagraph: you can only add a paragraph to a Section that ' ...
'contains only paragraphs']);
end
o.elements{end+1} = paragraph(varargin{:});
end

4
matlab/reports/@section/addTable.m
View File

@ -30,5 +30,9 @@ function o = addTable(o, varargin)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
for i=1:length(o.elements)
assert(~isa(o.elements{i}, 'paragraph'), ...
'@addTable: A Section that contains a Paratable cannot contain a Table');
end
o.elements{end+1} = report_table(varargin{:});
end

4
matlab/reports/@section/addVspace.m
View File

@ -29,5 +29,9 @@ function o = addVspace(o, varargin)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
for i=1:length(o.elements)
assert(~isa(o.elements{i}, 'paragraph'), ...
'@addVspace: A Section that contains a Paragraph cannot contain a Vspace');
end
o.elements{end+1} = vspace(varargin{:});
end

2
matlab/reports/@section/section.m
View File

@ -40,7 +40,7 @@ elseif nargin > 1
% overwrite default values
for pair = reshape(varargin, 2, [])
ind = strmatch(lower(pair{1}), lower(optNames), 'exact');
ind = find(strcmpi(optNames, pair{1}));
assert(isempty(ind) || length(ind) == 1);
if ~isempty(ind)
o.(optNames{ind}) = pair{2};

25
matlab/reports/@section/write.m
View File

@ -33,6 +33,7 @@ function o = write(o, fid, pg, sec)
assert(fid ~= -1);
fprintf(fid, '%% Section Object\n');
if ~isempty(o.height)
fprintf(fid, '\\setlength\\sectionheight{%s}%%\n', o.height);
end
@ -45,22 +46,34 @@ end
fprintf(fid, '}{%%\n');
fprintf(fid, '\\begin{tabular}[t]{');
for i=1:o.cols
fprintf(fid, 'c');
if ~isa(o.elements{1}, 'paragraph')
% if one element is a paragraph, they all are
% due to check in add*.m functions
fprintf(fid, 'l');
else
fprintf(fid, 'c');
end
end
fprintf(fid, '}\n');
ne = numElements(o);
nlcounter = 0;
row = 1;
col = 1;
for i=1:ne
if isa(o.elements{i}, 'vspace')
assert(col == o.cols, ['@section.write: must place ' ...
'vspace command after a linebreak in the table ' ...
'or series of charts']);
o.elements{i}.write(fid);
fprintf(fid, '\\\\\n');
if col ~= o.cols
fprintf(fid, '\\end{tabular}}\\\\\n');
fprintf(fid, '%% End Section Object\n\n');
return;
end
else
o.elements{i}.write(fid, pg, sec, row, col);
if isa(o.elements{i}, 'paragraph')
o.elements{i}.write(fid);
else
o.elements{i}.write(fid, pg, sec, row, col);
end
if col ~= o.cols
fprintf(fid, ' & ');
col = col + 1;

4
matlab/reports/@vspace/vspace.m
View File

@ -12,7 +12,7 @@ function o = vspace(varargin)
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2013 Dynare Team
% Copyright (C) 2013-2014 Dynare Team
%
% This file is part of Dynare.
%
@ -48,7 +48,7 @@ elseif nargin > 1
% overwrite default values
for pair = reshape(varargin, 2, [])
ind = strmatch(lower(pair{1}), lower(optNames), 'exact');
ind = find(strcmpi(optNames, pair{1}));
assert(isempty(ind) || length(ind) == 1);
if ~isempty(ind)
o.(optNames{ind}) = pair{2};

2
matlab/reports/allCellsAreDates.m
View File

@ -31,7 +31,7 @@ function tf = allCellsAreDates(dcell)
assert(iscell(dcell));
tf = true;
for i=1:length(dcell)
if ~isa(dcell{i}, 'dates')
if ~isdates(dcell{i})
tf = false;
return;
end

2
matlab/reports/allCellsAreDatesRange.m
View File

@ -31,7 +31,7 @@ function tf = allCellsAreDatesRange(dcell)
assert(iscell(dcell));
tf = true;
for i=1:length(dcell)
if ~(isa(dcell{i}, 'dates') && dcell{i}.ndat >= 2)
if ~(isdates(dcell{i}) && dcell{i}.ndat >= 2)
tf = false;
return;
end

10
matlab/reports/display_reporting_object.m
View File

@ -69,10 +69,14 @@ for i=1:length(fields)
fprintf('false');
end
elseif isobject(val)
if isa(val, 'dates')
fprintf('<dates: %s, ..., %s>', ...
if isdates(val)
if isempty(val)
fprintf('<dates: empty>');
else
fprintf('<dates: %s, ..., %s>', ...
date2string(val(1)), date2string(val(end)));
elseif isa(val, 'dseries')
end
elseif isdseries(val)
if numel(val) == 1
fprintf('<dseries: %s>', val.name{1});
else

14
matlab/reports/getMaxRange.m
View File

@ -23,12 +23,14 @@ ddmax = dates();
ne = length(cellser);
for i=1:ne
ddt = cellser{i}.getRange();
if isempty(ddmin)
ddmin = ddt(1);
ddmax = ddt(end);
else
ddmin = min(ddt(1), ddmin);
ddmax = max(ddt(end), ddmax);
if ~isempty(ddt)
if isempty(ddmin)
ddmin = ddt(1);
ddmax = ddt(end);
else
ddmin = min(ddt(1), ddmin);
ddmax = max(ddt(end), ddmax);
end
end
end
dd = ddmin:ddmax;

163
tests/ep/rs2.mod Normal file
View File

@ -0,0 +1,163 @@
var V lL Vkp Valphaexp lzn lzd lp0 lpi MC lwreal lY lDisp lC lpiavg Int lA lDZ lG pistar Intr
lsdf Int1 lpi1;
// pricebond pricebondrn ytm ytmrn termprem ehpr slope;
varexo epsA epsZ epsG epsPistar epsInt;
parameters theta xi beta phi alpha eta KBar chi0 LMax chi IBar rhoinflavg taylrho DZBar taylpi tayly YBar rhoa rhoz rhog rhopistar piBar gssload consoldelta VAIMSS GBar;
DZBar = 1.0025;
eta = 2/3;
IES = .11;
phi = 1/IES;
beta = .99 *DZBar^phi;
delta = .02;
Frisch = .28;
CRRA = 110;
LMax = 3;
// chi0 = exp(wrealAIMSS -phi*CAIMSS) *(LMax-1)^chi; // normalize L in ss = 1
chi0 = 1/3;
chi = 1/Frisch *(LMax-1);
alpha = (CRRA - 1/(1/phi + 1/chi *(LMax-1))) *(1/(1-phi) + 1/(1-chi) *(LMax-1));
theta = .2;
xi = .78;
K_Y = 10;
YBar = K_Y^((1-eta)/eta);
GBar = .17 *YBar;
KBar = K_Y *YBar;
IBar = (delta +(DZBar-1)) *KBar;
piBar = 1.0;
taylrho = .73;
taylpi = .53;
tayly = .93;
rhoa = .96;
rhoz = 0;
rhog = .95;
rhopistar = 0;
epsafl = 1; // these flags turn on or off each corresponding shock in the model
epsapermfl = 0;
epsgfl = 1;
epsifl = 1;
epspistarfl = 0;
rhoinflavg = .7; // for price level targeting; set rhoinflavg = .99 or .999
gssload = 0; // this is the GSS theta parameter in the long-run nominal risks section of the paper
VAIMSS = 1;
consoldelta = 1;
model;//(use_dll);
// Value function and Euler equation
V = exp(lC)^(1-phi) /(1-phi) + chi0 *(LMax-exp(lL))^(1-chi) /(1-chi) + beta *Vkp;
Int1 = Int/10;
lpi1 = lpi/10;
exp(lC)^-phi = beta *exp(10*Int1-10*lpi1(+1)) *exp(lC(+1))^-phi *exp(lDZ(+1))^-phi
*(V(+1) *exp(lDZ(+1))^(1-phi) /Vkp)^-alpha;
// The following two equations define the E-Z-W-K-P certainty equivalent term
// Vkp = (E_t V(+1)^(1-alpha))^(1/(1-alpha)). It takes two equations to do this because
// Perturbation AIM sets the expected value of all equations equal to zero; E_t F(variables) = 0.
// Thus; the first equation below defines Valphaexp = E_t V(+1)^(1-alpha). The second
// equation then takes the (1-alpha)th root of this expectation.
// Finally; the scaling and unscaling of Valphaexp by the constant VAIMSS and DZBar improves the
// numerical behavior of the model; without it; the steady-state value of Valphaexp can be minuscule
// (e.g.; 10^-50); which requires Mathematica to use astronomical levels of precision to solve. *)
[mcp = 'Valphaexp > 1e-5']
Valphaexp = (V(+1) *exp(lDZ(+1))^(1-phi) /VAIMSS /DZBar^(1-phi))^(1-alpha);
Vkp = VAIMSS *DZBar^(1-phi) *Valphaexp^(1/(1-alpha));
exp(lzn) = (1+theta) *MC *exp(lY) + xi *beta *exp(lC(+1)-lC)^-phi *exp(lDZ(+1))^-phi
*(V(+1) *exp(lDZ(+1))^(1-phi) /Vkp)^-alpha *exp(lpi(+1))^((1+theta)/theta/eta) *exp(lzn(+1));
exp(lzd) = exp(lY) + xi *beta *exp(lC(+1)-lC)^-phi*exp(lDZ(+1))^-phi
*(V(+1) *exp(lDZ(+1))^(1-phi) /Vkp)^-alpha *exp(lpi(+1))^(1/theta) *exp(lzd(+1));
exp(lp0)^(1+(1+theta)/theta *(1-eta)/eta) = exp(lzn-lzd);
exp(lpi)^(-1/theta) = (1-xi) *exp(lp0+lpi)^(-1/theta) + xi;
// Marginal cost and real wage
MC = exp(lwreal) /eta *exp(lY)^((1-eta)/eta) /exp(lA)^(1/eta) /KBar^((1-eta)/eta);
[mcp = 'lL < 1.0986']
chi0 *(LMax-exp(lL))^-chi /exp(lC)^-phi = exp(lwreal);
// Output equations
exp(lY) = exp(lA) *KBar^(1-eta) *exp(lL)^eta /exp(lDisp);
exp(lDisp)^(1/eta) = (1-xi) *exp(lp0)^(-(1+theta)/theta/eta)
+ xi *exp(lpi)^((1+theta)/theta/eta) *exp(lDisp(-1))^(1/eta);
exp(lC) = exp(lY)-exp(lG)-IBar; // aggregate resource constraint; no adj costs
// Monetary Policy Rule
lpiavg = rhoinflavg *lpiavg(-1) + (1-rhoinflavg) *lpi;
4*Int = (1-taylrho) * ( 4*log(1/beta *DZBar^phi) + 4*lpiavg
+ taylpi * (4*lpiavg-pistar) + tayly * (exp(lY)-YBar)/YBar )
+ taylrho * 4*Int(-1) + epsInt; // multiply Int; infl by 4 to put at annual rate
// Exogenous Shocks
lA = rhoa * lA(-1) + epsA;
lDZ = (1-rhoz)*log(DZBar) + rhoz * lDZ(-1) + epsZ;
lG = (1-rhog)*log(GBar) + rhog * lG(-1) + epsG;
pistar = (1-rhopistar) *log(piBar) + rhopistar *pistar(-1) + gssload *(4*lpiavg-pistar)
+ epsPistar;
// Term premium and other auxiliary finance equations
Intr = Int(-1) - lpi; // ex post real short rate
exp(lsdf) = beta *exp(lC(+1)-lC)^-phi *exp(lDZ(+1))^-phi
*(V(+1) *exp(lDZ(+1))^(1-phi) /Vkp)^-alpha /exp(lpi(+1));
//pricebond = 1 + consoldelta *beta *exp(lC(+1)-lC)^-phi *exp(lDZ(+1))^-phi
// *(V(+1) *exp(lDZ(+1))^(1-phi) /Vkp)^-alpha /pi(+1) *pricebond(+1);
//pricebondrn = 1 + consoldelta *pricebondrn(+1) /exp(Int);
//ytm = log(consoldelta*pricebond/(pricebond-1)) *400; // yield in annualized pct
//ytmrn = log(consoldelta*pricebondrn/(pricebondrn-1)) *400;
//termprem = 100 * (ytm–ytmrn); // term prem in annualized basis points
//ehpr = ( (consoldelta *pricebond + exp(Int(-1))) /pricebond(-1)–exp(Int(-1)) *400;
//slope = ytm–Int*400;
end;
steady_state_model;
lA = 0;
lDZ = log(DZBar);
lG = log(GBar);
pistar = log(piBar);
lpi = log(piBar);
lpiavg = lpi;
lDisp = 0;
Int = log(DZBar^phi/beta)+lpi;
Int1 = Int/10;
lpi1 = lpi/10;
lL = 0;
lY = log(YBar);
lC = log(exp(lY)-exp(lG)-IBar);
lp0 = log(((exp(lpi)^(-1/theta)-xi)/((1-xi)*exp(lpi)^(-1/theta)))^(-theta));
lzd = log(exp(lY)/(1-xi *beta *exp(lDZ)^-phi
*exp(lpi)^(1/theta)));
lzn = log(exp(lp0)^(1+(1+theta)/theta *(1-eta)/eta)*exp(lzd));
MC = (exp(lzn)-xi *beta *DZBar^-phi*piBar^((1+theta)/theta/eta)*exp(lzn))/((1+theta)*YBar);
lwreal = log(MC*eta /(exp(lY)^((1-eta)/eta) / KBar^((1-eta)/eta)));
chi0 = exp(lwreal)*exp(lC)^(-phi)/(LMax-exp(lL))^(-chi);
V = (exp(lC)^(1-phi) /(1-phi) + chi0 *(LMax-exp(lL))^(1-chi) /(1-chi))/(1-beta*DZBar^(1-phi));
k = 1.05;
VAIMSS = k*V;
Valphaexp = k^(alpha-1);
Vkp = DZBar^(1-phi)*V;
Intr = Int - lpi;
lsdf = log(beta *exp(lDZ)^-phi/exp(lpi));
end;
steady;
shocks;
var epsA; stderr 0.005;
var epsZ; stderr 0.001;
var epsG; stderr 0.004;
var epsPistar; stderr 0.0005;
var epsInt; stderr 0.003;
end;
//stoch_simul(order=3,periods=50000,pruning);
options_.ep.verbosity=0;
options_.ep.stochastic.algo=1;
options_.ep.solve_algo = 10;
options_.ep.maxit = 100;
options_.ep.IntegrationAlgorithm='UT_2p+1';
options_.ep.ut.k = 1;
options_.solve_tolf = 1e-12;
extended_path(order=1,periods=3);

31
tests/reporting/runDynareReport.m
View File

@ -144,7 +144,7 @@ rep = rep.addSection();
rep = CommResidTablePage(rep, db_q, dc_q, trange, dates('2012q2'));
%% Commodities Graphs
%Page 1
%Page 24
rep = rep.addPage('title', 'Jan1 vs Jan2', ...
'titleFormat', '\large\bfseries');
rep = rep.addSection();
@ -184,12 +184,25 @@ rep = rep.addSeries('data', db_q{'LRPFOOD_BAR_WORLD'}, ...
'graphLineColor', 'green', ...
'graphLineStyle', 'solid', ...
'graphLineWidth', 1.5);
rep = rep.addSeries('graphVline', dates('2009q2'), ...
'graphLineColor', 'red', ...
'graphLineWidth', 1.5);
% Pae 2
% Page 25
rep = rep.addPage('title', {'Jan1 vs Jan2', 'World Oil and Food Prices'}, ...
'titleFormat', {'\large\bfseries', '\large'});
rep = rep.addSection('cols', 2);
rep = rep.addSection('cols', 1);
rep = rep.addParagraph('text', 'Lorem ipsum dolor sit amet, consectetur adipisicing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua. Ut enim ad minim veniam, quis nostrud exercitation ullamco laboris nisi ut aliquip ex ea commodo consequat. Duis aute irure dolor in reprehenderit in voluptate velit esse cillum dolore eu fugiat nulla pariatur. Excepteur sint occaecat cupidatat non proident, sunt in culpa qui officia deserunt mollit anim id est laborum.', ...
'cols', 2, ...
'heading', '\textbf{My First Paragraph Has Two Columns}');
rep = rep.addSection('cols', 1);
rep = rep.addParagraph('text', 'Lorem ipsum dolor sit amet, consectetur adipisicing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua. Ut enim ad minim veniam, quis nostrud exercitation ullamco laboris nisi ut aliquip ex ea commodo consequat.\newline', ...
'heading', '\textbf{My Next Paragraphs Only Have One}', ...
'indent', false);
rep = rep.addParagraph('text', 'Lorem ipsum dolor sit amet, consectetur adipisicing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua. Ut enim ad minim veniam, quis nostrud exercitation ullamco laboris nisi ut aliquip ex ea commodo consequat. Duis aute irure dolor in reprehenderit in voluptate velit esse cillum dolore eu fugiat nulla pariatur. Excepteur sint occaecat cupidatat non proident, sunt in culpa qui officia deserunt mollit anim id est laborum.\newline');
rep = rep.addSection('cols', 2);
rep = rep.addGraph('title', 'World Real Oil Price', ...
'xrange', prange, ...
@ -231,6 +244,10 @@ rep = rep.addSeries('data', dc_q{'LRPFOOD_WORLD'}, ...
'graphLineColor', 'blue', ...
'graphLineStyle', 'dashed', ...
'graphLineWidth', 1.5);
rep = rep.addSeries('graphHline', 460, ...
'graphLineColor', 'red', ...
'graphLineWidth', 1.5);
rep = rep.addGraph('title', 'Equilibrium World Real Food Price', ...
'xrange', prange, ...
@ -239,11 +256,15 @@ rep = rep.addGraph('title', 'Equilibrium World Real Food Price', ...
'xTickLabelRotation', 0);
rep = rep.addSeries('data', db_q{'LRPFOOD_BAR_WORLD'}, ...
'graphLineColor', 'blue', ...
'graphLineWidth', 1.5);
'graphLineWidth', 1.5, ...
'graphLegendName', 'baseline', ...
'graphMiscTikzAddPlotOptions', 'mark=halfcircle*,color=red');
rep = rep.addSeries('data', dc_q{'LRPFOOD_BAR_WORLD'}, ...
'graphLineColor', 'blue', ...
'graphLineStyle', 'dashed', ...
'graphLineWidth', 1.5);
'graphLineWidth', 1.5, ...
'graphLegendName', 'control', ...
'graphMiscTikzAddPlotOptions', 'mark=halfcircle*,mark options={rotate=90,scale=3}');
%% Write & Compile Report
rep.write();