Add example for system prior

(@stepan removed trailing spaces)

(cherry picked from commit cf6621187e)

.gitlab-ci.yml

@@ -10,12 +10,10 @@ variables:
 # - if VERSION was already set (when manually running a pipeline), use it
 # - if we are in the official Dynare repository:
 #   + if on a tag: use the tag
-#   + if on master: use 4.7-unstable-$TIMESTAMP-$COMMIT
-#   + on another branch: use $BRANCH-$TIMESTAMP-$COMMIT
+#   + otherwise: use $BRANCH-$TIMESTAMP-$COMMIT
 # - if in a personal repository: use $USER-$TIMESTAMP-$COMMIT
 before_script:
   - '[[ -z $VERSION ]] && [[ $CI_PROJECT_NAMESPACE == Dynare ]] && [[ -n $CI_COMMIT_TAG ]] && export VERSION=$CI_COMMIT_TAG'
-  - '[[ -z $VERSION ]] && [[ $CI_PROJECT_NAMESPACE == Dynare ]] && [[ $CI_COMMIT_REF_NAME == master ]] && export VERSION=4.7-unstable-$(date +%F-%H%M)-$CI_COMMIT_SHORT_SHA'
   - '[[ -z $VERSION ]] && [[ $CI_PROJECT_NAMESPACE == Dynare ]] && export VERSION=$CI_COMMIT_REF_NAME-$(date +%F-%H%M)-$CI_COMMIT_SHORT_SHA'
   - '[[ -z $VERSION ]] && export VERSION=$CI_PROJECT_NAMESPACE-$(date +%F-%H%M)-$CI_COMMIT_SHORT_SHA'
 
@@ -184,10 +182,10 @@ test_dynare++:
 # those jobs to start before the “test” and “pkg” stages have succeeded. Hence
 # we stick to the “dependencies” keyword.
 
-deploy_manual_unstable:
+deploy_manual_stable:
   stage: deploy
   rules:
-    - if: '$CI_PROJECT_NAMESPACE == "Dynare" && $CI_COMMIT_REF_NAME == "master"'
+    - if: '$CI_PROJECT_NAMESPACE == "Dynare" && $CI_COMMIT_TAG =~ /^4\.7\.[0-9]+$/'
       when: on_success
     - when: never
   tags:
@@ -197,12 +195,13 @@ deploy_manual_unstable:
   script:
     - rm -rf doc/manual/build/html/_static/mathjax
     - ln -s /usr/share/javascript/mathjax doc/manual/build/html/_static/mathjax
-    - rsync --recursive --links --delete doc/manual/build/html/ /srv/www.dynare.org/manual-unstable/
+    - rsync --recursive --links --delete doc/manual/build/html/ /srv/www.dynare.org/manual/
+    - cp doc/manual/build/latex/dynare-manual.pdf /srv/www.dynare.org/manual.pdf
 
-deploy_snapshot_unstable:
+deploy_release_stable:
   stage: deploy
   rules:
-    - if: '$CI_PROJECT_NAMESPACE == "Dynare" && $CI_COMMIT_REF_NAME == "master"'
+    - if: '$CI_PROJECT_NAMESPACE == "Dynare" && $CI_COMMIT_TAG =~ /^4\.7\.[0-9]+$/'
       when: on_success
     - when: never
   tags:
@@ -213,10 +212,30 @@ deploy_snapshot_unstable:
     - pkg_macOS
   script:
     - f=(windows/exe/*) && osslsigncode sign -pkcs12 ~/cepremap-comodo-sectigo-code-signing.p12 -n Dynare -i https://www.dynare.org -t http://timestamp.comodoca.com -in ${f[0]} -out ${f[0]}.signed && mv ${f[0]}.signed ${f[0]}
-    - cp *.tar.xz /srv/www.dynare.org/snapshot/source/ && ln -sf *.tar.xz /srv/www.dynare.org/snapshot/source/dynare-latest-src.tar.xz
-    - f=(windows/exe/*) && cp ${f[0]} /srv/www.dynare.org/snapshot/windows/ && ln -sf ${f[0]##*/} /srv/www.dynare.org/snapshot/windows/dynare-latest-win.exe
-    - f=(windows/7z/*) && cp ${f[0]} /srv/www.dynare.org/snapshot/windows-7z/ && ln -sf ${f[0]##*/} /srv/www.dynare.org/snapshot/windows-7z/dynare-latest-win.7z
-    - f=(windows/zip/*) && cp ${f[0]} /srv/www.dynare.org/snapshot/windows-zip/ && ln -sf ${f[0]##*/} /srv/www.dynare.org/snapshot/windows-zip/dynare-latest-win.zip
-    - f=(macOS/pkg/*) && cp ${f[0]} /srv/www.dynare.org/snapshot/macos/ && ln -sf ${f[0]##*/} /srv/www.dynare.org/snapshot/macos/dynare-latest-macos.pkg
-    - ~/update-snapshot-list.sh
+    - cp *.tar.xz /srv/www.dynare.org/release/source/
+    - cp windows/exe/* /srv/www.dynare.org/release/windows/
+    - cp windows/7z/* /srv/www.dynare.org/release/windows-7z/
+    - cp windows/zip/* /srv/www.dynare.org/release/windows-zip/
+    - cp macOS/pkg/* /srv/www.dynare.org/release/macos/
+    - ~/update-release-list.sh
     - curl -X POST -F token="$WEBSITE_PIPELINE_TRIGGER_TOKEN" -F ref=master https://git.dynare.org/api/v4/projects/40/trigger/pipeline
+
+deploy_beta_stable:
+  stage: deploy
+  rules:
+    - if: '$CI_PROJECT_NAMESPACE == "Dynare" && $CI_COMMIT_TAG =~ /^4\.7(\.[0-9]+)?-(beta|rc)[0-9]+$/'
+      when: on_success
+    - when: never
+  tags:
+    - restricted
+  dependencies:
+    - pkg_source
+    - pkg_windows
+    - pkg_macOS
+  script:
+    - f=(windows/exe/*) && osslsigncode sign -pkcs12 ~/dynare-object-signing.p12 -n Dynare -i https://www.dynare.org -t http://timestamp.digicert.com -in ${f[0]} -out ${f[0]}.signed && mv ${f[0]}.signed ${f[0]}
+    - cp *.tar.xz /srv/www.dynare.org/beta/source/
+    - cp windows/exe/* /srv/www.dynare.org/beta/windows/
+    - cp windows/7z/* /srv/www.dynare.org/beta/windows-7z/
+    - cp windows/zip/* /srv/www.dynare.org/beta/windows-zip/
+    - cp macOS/pkg/* /srv/www.dynare.org/beta/macos/

configure.ac

@@ -1,6 +1,6 @@
 dnl Process this file with autoconf to produce a configure script.
 
-dnl Copyright © 2009-2019 Dynare Team
+dnl Copyright © 2009-2021 Dynare Team
 dnl
 dnl This file is part of Dynare.
 dnl
@@ -18,7 +18,7 @@ dnl You should have received a copy of the GNU General Public License
 dnl along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
 AC_PREREQ([2.62])
-AC_INIT([dynare], [4.7-unstable])
+AC_INIT([dynare], [4.7.0])
 AC_CONFIG_SRCDIR([preprocessor/src/DynareMain.cc])
 AM_INIT_AUTOMAKE([1.11 -Wall -Wno-portability foreign no-dist-gzip dist-xz tar-pax])
 

doc/bvar-a-la-sims.tex

@@ -21,7 +21,7 @@
     at: \url{https://www.gnu.org/licenses/fdl.txt}
     \newline
     \indent Many thanks to Christopher Sims for providing his BVAR
-    MATLAB\textregistered~routines, to St\'ephane Adjemian and Michel Juillard
+    MATLAB\textsuperscript{\textregistered}~routines, to St\'ephane Adjemian and Michel Juillard
     for their helpful support, and to Marek Jaroci\'nski for reporting a bug.
   }}
 

doc/dr.bib

@@ -1,111 +1,102 @@
+% Encoding: UTF-8
 
-@techreport{adjemian/al:2011,
-	author = {Adjemian, St\'ephane and Bastani, Houtan and Juillard, Michel and Mihoubi, Ferhat and Perendia, George and Ratto, Marco and Villemot, S\'ebastien},
-	title = {Dynare: Reference Manual, Version 4},
-	institution = {CEPREMAP},
-	year = {2011},
-	type = {Dynare Working Papers},
-	number = {1}
+@TechReport{adjemian/al:2011,
+  author      = {Adjemian, St\'ephane and Bastani, Houtan and Juillard, Michel and Karam\'e, Fr\'ederic and Maih, Junior and Mihoubi, Ferhat and Mutschler, Willi and Perendia, George and Pfeifer, Johannes and Ratto, Marco and Villemot, S\'ebastien},
+  institution = {CEPREMAP},
+  title       = {Dynare: Reference Manual Version 4},
+  year        = {2011},
+  number      = {1},
+  type        = {Dynare Working Papers},
+}
+
+@Article{blanchard/kahn:1980,
+  author   = {Blanchard, Olivier Jean and Kahn, Charles M.},
+  journal  = {Econometrica},
+  title    = {The Solution of Linear Difference Models under Rational Expectations},
+  year     = {1980},
+  month    = {7},
+  number   = {5},
+  pages    = {1305-11},
+  volume   = {48},
+  abstract = {Many have questioned the empirical relevance of the Calvo-Yun model. This paper adds a term structure to three widely studied macroeconomic models (Calvo-Yun, hybrid and Svensson). We back out from observations on the yield curve the underlying macroeconomic model that most closely matches the level, slope and curvature of the yield curve. With each model we trace the response of the yield curve to macroeconomic shocks. We assess the fit of each model against the observed behaviour of interest rates and find limited support for the Calvo-Yun model in terms of fit with the observed yield curve, we find some support for the hybrid model but the Svensson model performs best.},
+  doi      = {10.2307/1912186},
+  keywords = {Macromodels Yield curve Persistence},
+}
+
+@Article{klein:2000,
+  author  = {Klein, Paul},
+  journal = {Journal of Economic Dynamics and Control},
+  title   = {Using the generalized {Schur} form to solve a multivariate linear rational expectations model},
+  year    = {2000},
+  month   = {September},
+  number  = {10},
+  pages   = {1405-1423},
+  volume  = {24},
+  doi     = {10.1016/s0165-1889(99)00045-7},
+}
+
+@Article{schmitt-grohe/uribe:2004,
+  author  = {Schmitt-Groh\'{e}, Stephanie and Ur\'{i}be, Martin},
+  journal = {Journal of Economic Dynamics and Control},
+  title   = {Solving dynamic general equilibrium models using a second-order approximation to the policy function},
+  year    = {2004},
+  month   = {January},
+  number  = {4},
+  pages   = {755-775},
+  volume  = {28},
+  doi     = {10.1016/s0165-1889(03)00043-5},
+}
+
+@Article{sims:2001,
+  author  = {Sims, Christopher A},
+  journal = {Computational Economics},
+  title   = {Solving Linear Rational Expectations Models},
+  year    = {2002},
+  month   = {October},
+  number  = {1-2},
+  pages   = {1-20},
+  volume  = {20},
+  doi     = {10.1023/A:1020517101123},
+}
+
+@InCollection{uhlig:1999,
+  author    = {Uhlig, Harald},
+  booktitle = {Computational Methods for the Study of Dynamic Economies},
+  publisher = {Oxford University Press},
+  address   = {Oxford},
+  title     = {A toolkit for analysing nonlinear dynamic stochastic models easily},
+  year      = {1999},
+  editor    = {Marimon, Ramon and Scott, Andrew},
+  pages     = {30-61},
 }
 
 
-@article{blanchard/kahn:1980,
-	author = {Blanchard, Olivier Jean and Kahn, Charles M.},
-	title = {The Solution of Linear Difference Models under Rational Expectations},
-	journal = {Econometrica},
-	year = 1980,
-	volume = {48},
-	number = {5},
-	pages = {1305-11},
-	month = {July},
-	keywords = { Macromodels Yield curve Persistence},
-	abstract = {Many have questioned the empirical relevance of the Calvo-Yun model. This paper adds a term structure to three widely studied macroeconomic models (Calvo-Yun, hybrid and Svensson). We back out from observations on the yield curve the underlying macroeconomic model that most closely matches the level, slope and curvature of the yield curve. With each model we trace the response of the yield curve to macroeconomic shocks. We assess the fit of each model against the observed behaviour of interest rates and find limited support for the Calvo-Yun model in terms of fit with the observed yield curve, we find some support for the hybrid model but the Svensson model performs best.},
-	url = {http://ideas.repec.org/a/ecm/emetrp/v48y1980i5p1305-11.html}
-}
-
-
-@article{klein:2000,
-	author = {Klein, Paul},
-	title = {Using the generalized Schur form to solve a multivariate linear rational expectations model},
-	journal = {Journal of Economic Dynamics and Control},
-	year = 2000,
-	volume = {24},
-	number = {10},
-	pages = {1405-1423},
-	month = {September},
-	keywords = {},
-	abstract = {},
-	url = {http://ideas.repec.org/a/eee/dyncon/v24y2000i10p1405-1423.html}
-}
-
-
-@article{schmitt-grohe/uribe:2004,
-	author = {Schmitt-Groh\'{e}, Stephanie and Ur\'{i}be, Martin},
-	title = {Solving dynamic general equilibrium models using a second-order approximation to the policy function},
-	journal = {Journal of Economic Dynamics and Control},
-	year = 2004,
-	volume = {28},
-	number = {4},
-	pages = {755-775},
-	month = {January},
-	keywords = {},
-	url = {http://ideas.repec.org/a/eee/dyncon/v28y2004i4p755-775.html}
-}
-
-
-@article{sims:2001,
-	author = {Sims, Christopher A},
-	title = {Solving Linear Rational Expectations Models},
-	journal = {Computational Economics},
-	year = 2002,
-	volume = {20},
-	number = {1-2},
-	pages = {1-20},
-	month = {October},
-	keywords = {},
-	abstract = {},
-	url = {http://ideas.repec.org/a/kap/compec/v20y2002i1-2p1-20.html}
-}
-
-
-@incollection{uhlig:1999,
-	author = {Uhlig, Harald},
-	title = {A toolkit for analysing nonlinear dynamic stochastic models easily},
-	booktitle = {Computational Methods for the Study of Dynamic Economics},
-	publisher = {Oxford University Press},
-	year = {1999},
-	editor = {Marimon, Ramon and Scott, Androw},
-	pages = {30-61}
-}
-
-
-@techreport{kamenik:2003,
+@techreport{kamenik:2004,
 	author = {Kamenik, Ondra},
 	title = {Solution of Specialized Sylvester Equation},
 	institution = {Manuscript},
-	year = {2003}
+	year = {2004}
 }
 
-
-@article{collard/juillard:2001:compecon,
-	author = {Collard, Fabrice and Juillard, Michel},
-	title = {A Higher-Order Taylor Expansion Approach to Simulation of Stochastic Forward-Looking Models with an Application to a Nonlinear Phillips Curve Model},
-	journal = {Computational Economics},
-	year = {2001},
-	volume = {17},
-	number = {2-3},
-	pages = {125-39},
-	month = {June},
-	keywords = {},
-	url = {http://ideas.repec.org/a/kap/compec/v17y2001i2-3p125-39.html}
+@Article{collard/juillard:2001:compecon,
+  author  = {Collard, Fabrice and Juillard, Michel},
+  journal = {Computational Economics},
+  title   = {A Higher-Order {Taylor} Expansion Approach to Simulation of Stochastic Forward-Looking Models with an Application to a Nonlinear {Phillips} Curve Model},
+  year    = {2001},
+  month   = {6},
+  number  = {2-3},
+  pages   = {125-139},
+  volume  = {17},
+  doi     = {10.1023/A:1011624124377},
 }
 
-
-@book{golub/van-loan:1996,
-	author = {Golub, Gene H. and Van Loan, Charles F.},
-	title = {Matrix Computations},
-	publisher = {The John Hopkins University Press},
-	year = {1996},
-	edition = {third}
+@Book{golub/van-loan:1996,
+  author    = {Golub, Gene H. and Van Loan, Charles F.},
+  publisher = {The John Hopkins University Press},
+  title     = {Matrix Computations},
+  year      = {2013},
+  address   = {Baltimore},
+  edition   = {4},
 }
 
+@Comment{jabref-meta: databaseType:bibtex;}

doc/dr.tex

@@ -3,7 +3,7 @@
 \usepackage{amssymb}
 \usepackage{amsmath}
 \usepackage{hyperref}
-\hypersetup{breaklinks=true,pagecolor=white,colorlinks=true,linkcolor=blue,citecolor=blue,urlcolor=blue}
+\hypersetup{breaklinks=true,colorlinks=true,linkcolor=blue,citecolor=blue,urlcolor=blue}
 \usepackage{natbib}
 
 \usepackage{fullpage}
@@ -36,7 +36,7 @@
   computing the first order approximated solution of a nonlinear rational
   expectations model. The core of the algorithm is a generalized Schur
   decomposition (also known as the QZ decomposition), as advocated by several
-  authors in the litterature. The contribution of the present paper is to focus
+  authors in the literature. The contribution of the present paper is to focus
   on implementation details that make the algorithm more generic and more
   efficient, especially for large models.
 
@@ -62,12 +62,8 @@ detail the algorithm implemented in Dynare for computing the first order
 approximated solution of nonlinear rational expectations models.\footnote{This
   algorithm is available using the \texttt{stoch\_simul} command of Dynare. The
   original implementation of this algorithm was done by Michel Juillard, using
-  MATLAB, and is available in the \texttt{matlab/dr1.m} file which is
-  distributed with Dynare. Another implementation was done by the author, in
-  C++, in the \texttt{DecisionRules} class, in the
-  \texttt{mex/sources/estimation} directory of the source tree. The notations
-  used in the present paper are closer to the C++ implementation than to the
-  MATLAB implementation.}
+  MATLAB, and is available in the \texttt{matlab/dyn\_first\_order\_solver.m} file which is
+  distributed with Dynare.}
 
 This algorithm is based on a generalized Schur decomposition---also known as
 the QZ decomposition---and is therefore essentially a variation on the
@@ -653,8 +649,8 @@ equation is:
 \end{equation*}
 In the general case, this equation is a specialized Sylvester equation, which
 can be solved using the algorithm proposed by
-\citet{kamenik:2003}\footnote{This paper is distributed with Dynare, in the
-  \texttt{sylvester.pdf} file under the documentation directory.}.
+\citet{kamenik:2004}\footnote{This paper is distributed with Dynare, in the
+  \texttt{sylvester.pdf} file under the Dynare++ documentation directory.}.
 
 \bibliographystyle{elsarticle-harv}
 \bibliography{dr}

doc/gsa/gsa.tex

@@ -4,15 +4,15 @@
 \documentclass[12pt,a4paper]{article}
 \usepackage{amssymb,amsmath}
 \usepackage[dvips]{graphicx}
-\usepackage{natbib}
 \usepackage{psfrag}
 \usepackage{setspace}
 \usepackage{rotating}
 \usepackage{hyperref}
-\hypersetup{breaklinks=true,pagecolor=white,colorlinks=true,linkcolor=blue,citecolor=blue,urlcolor=blue}
+\hypersetup{breaklinks=true,colorlinks=true,linkcolor=blue,citecolor=blue,urlcolor=blue}
 %\singlespacing (interlinea singola)
 %\onehalfspacing (interlinea 1,5)
 %\doublespacing (interlinea doppia)
+\usepackage{doi,natbib}
 
 
 %\bibpunct{(}{)}{;}{a}{,}{,}

doc/gsa/marco.bib

@@ -1,19 +1,24 @@
-@ARTICLE{Ratto_CompEcon_2008,
-  author =       {Ratto, M.},
-  title =        {Analysing DSGE Models with Global Sensitivity Analysis},
-  journal =      {Computational Economics},
-  year =         {2008},
-  volume =       {31},
-  pages =        {115--139},
+% Encoding: UTF-8
+
+@Article{Ratto_CompEcon_2008,
+  author  = {Ratto, Marco},
+  journal = {Computational Economics},
+  title   = {Analysing {DSGE} models with global sensitivity analysis},
+  year    = {2008},
+  pages   = {115--139},
+  volume  = {31},
+  doi     = {10.1007/s10614-007-9110-6},
 }
 
-@ARTICLE{Iskrev2010,
-  author = {Nikolay Iskrev},
-  title = {Local Identification in {DSGE} Models},
+@Article{Iskrev2010,
+  author  = {Nikolay Iskrev},
   journal = {Journal of Monetary Economics},
-  year = {2010},
-  volume = {57},
-  pages = {189-202}
+  title   = {Local Identification in {DSGE} Models},
+  year    = {2010},
+  number  = {2},
+  pages   = {189-202},
+  volume  = {57},
+  doi     = {10.1016/j.jmoneco.2009.12.007},
 }
 
 @UNPUBLISHED{Iskrev2011,
@@ -23,3 +28,5 @@
   note = {mimeo},
   year = {2011}
 }
+
+@Comment{jabref-meta: databaseType:bibtex;}

doc/guide.tex

@@ -1,23 +1,23 @@
 \documentclass[11pt,a4paper]{article}
-\usepackage{bibmad,graphicx,latexsym,amssymb,times}
-\usepackage[cp850]{inputenc}
+\usepackage{graphicx,latexsym,amssymb,times}
+\usepackage[utf8]{inputenc}
 
 \begin{document}
 \title{Stochastic simulations with {\sc Dynare}. \\ A practical guide.}
-\author{Fabrice Collard (GREMAQ, University of Toulouse)\\Adapted for Dynare 4.1\\ by Michel Juillard and S\'ebastien Villemot (CEPREMAP)}
+\author{Fabrice Collard (GREMAQ, University of Toulouse)\\Adapted for Dynare 4.x\\ by Michel Juillard and S\'ebastien Villemot (CEPREMAP)}
 \date{First draft: February 2001\hspace{10mm}This draft: December 2009.}
 \maketitle
 This document describes a model involving both endogenous and exogenous state variable. We first describe the theoretical model, before showing how the perturbation method is implemented in {\sc Dynare}.
 
-\section{A theoretical model} 
-We consider an economy that consists of a large number of dynastic households and a large number of firms. Firms are producing a homogeneous final product that can be either consumed or invested by means of capital and labor services. Firms own their capital stock and hire labor supplied by the households. Households own the firms. In each and every period three perfectly competitive markets open --- the markets for consumption goods, labor services, and financial capital in the form of firms' shares. 
+\section{A theoretical model}
+We consider an economy that consists of a large number of dynastic households and a large number of firms. Firms are producing a homogeneous final product that can be either consumed or invested by means of capital and labor services. Firms own their capital stock and hire labor supplied by the households. Households own the firms. In each and every period three perfectly competitive markets open --- the markets for consumption goods, labor services, and financial capital in the form of firms' shares.
 Household preferences are characterized by the lifetime utility function:
 \begin{equation}
   E_t\sum_{\tau=t}^{\infty}{\beta^\star}^{\tau-t} \left(\log(c_t)-\theta\frac{h_t^{1+\psi}}{1+\psi}\right)
   \label{eq:ut}
 \end{equation}
 \noindent where $0<\beta^\star<1$ is a constant discount factor, $c_t$ is consumption in period
-$t$, $h_t$ is the fraction of total available time devoted to productive activity in period $t$, $\theta>0$ and $\psi\geqslant 0$. We assume that there exists a central planner that determines hours, consumption and capital accumulation maximizing the household's utility function subject to the following budget constraint 
+$t$, $h_t$ is the fraction of total available time devoted to productive activity in period $t$, $\theta>0$ and $\psi\geqslant 0$. We assume that there exists a central planner that determines hours, consumption and capital accumulation maximizing the household's utility function subject to the following budget constraint
 \begin{equation}
   c_t+i_t=y_t
   \label{eq:bud}
@@ -44,7 +44,7 @@ a_t\\b_t
 \left(
 \begin{array}{cc}
 \rho&\tau\\
-\tau&\rho\\ 
+\tau&\rho\\
 \end{array}
 \right)\left(
 \begin{array}{c}
@@ -57,7 +57,7 @@ a_{t-1}\\b_{t-1}
 \end{array}
 \right)  \label{eq:process}
 \end{equation}
-where $|\rho+\tau|<1$ and $|\rho-\tau|<1 $ for sake of stationarity and 
+where $|\rho+\tau|<1$ and $|\rho-\tau|<1 $ for sake of stationarity and
 \begin{eqnarray*}
 E(\varepsilon_t)&=& 0,\\
 E(\nu_t)&=& 0,\\
@@ -65,7 +65,7 @@ E(\varepsilon_t\varepsilon_s)&=&\left\{
 \begin{array}{lcl}
 \sigma^2_\varepsilon & \mbox{ if } & t=s \\
 0 & \mbox{ if } & t\neq s \\
-\end{array}\right. \mbox{, }\\ 
+\end{array}\right. \mbox{, }\\
 E(\nu_t\nu_s)&=&\left\{
 \begin{array}{lcl}
 \sigma^2_\nu & \mbox{ if } & t=s \\
@@ -75,24 +75,24 @@ E(\varepsilon_t\nu_s)&=&\left\{
 \begin{array}{lcl}
 \varphi\sigma_\varepsilon\sigma_\nu & \mbox{ if } & t=s \\
 0 & \mbox{ if } & t\neq s \\
-\end{array}\right. \mbox{. } 
+\end{array}\right. \mbox{. }
 \end{eqnarray*}
 
- 
+
 \section{Dynamic Equilibrium}
-The dynamic equilibrium of this economy follows from the first order conditions for optimality: 
+The dynamic equilibrium of this economy follows from the first order conditions for optimality:
 \begin{eqnarray*}
 &&c_t \theta h_t^{1+\psi}=(1-\alpha) y_t \\
 &&\beta E_t\left[\left(\frac{\exp(b_t) c_t}{\exp(b_{t+1})c_{t+1}}\right)\left(\exp(b_{t+1})\alpha \frac{y_{t+1}}{k_{t+1}}+1-\delta\right)\right]=1\\
 &&y_t=\exp(a_t) k_t^\alpha h_t^{1-\alpha} \\
 &&k_{t+1}=\exp(b_t)(y_t-c_t)+(1-\delta)k_t \\
 &&a_t=\rho a_{t-1}+\tau b_{t-1}+\varepsilon_t \\
-&&b_t=\tau a_{t-1}+\rho b_{t-1}+\nu_t 
+&&b_t=\tau a_{t-1}+\rho b_{t-1}+\nu_t
 \end{eqnarray*}
 \section{The {\sc dynare} code}
 The {\sc dynare} code is straightforward to write, as the equilibrium is written in the natural way. The whole code is reported at the end of the section. Before that we proceed step by step.
 \paragraph{Preamble}
-The preamble consists of the some declarations to setup the endogenous and exogenous variables, the parameters and assign values to these parameters. 
+The preamble consists of the some declarations to setup the endogenous and exogenous variables, the parameters and assign values to these parameters.
 \begin{enumerate}
 \item {\tt var y, c, k, h, a, b;} specifies the endogenous variables in the model since we have output ({\tt y}), consumption ({\tt c}), capital ({\tt k}), hours ({\tt h}) and the two shocks ({\tt a, b}).
 \item {\tt varexo e, u;} specifies the exogenous variables in the model --- namely the innovations of the shocks, since we have the innovation of the non--incorporated shock ({\tt e}), and the innovation of the incorporated shock ({\tt u}).
@@ -133,9 +133,9 @@ theta = 2.95;
 \end{enumerate}
 
 \paragraph{Declaration of the model:}
-This step is done in a straightforward way. It starts with the instruction {\tt model;} and ends with {\tt end;}, in between all equilibrium conditions are written exactly the way we write it ``by hand''. However, there is a simple rule that should be kept in mind when the model is written. Let us consider a variable $x$: 
+This step is done in a straightforward way. It starts with the instruction {\tt model;} and ends with {\tt end;}, in between all equilibrium conditions are written exactly the way we write it ``by hand''. However, there is a simple rule that should be kept in mind when the model is written. Let us consider a variable $x$:
 \begin{itemize}
-\item If $x$ is decided in period $t$ then we simply write ${\tt x}$. 
+\item If $x$ is decided in period $t$ then we simply write ${\tt x}$.
 \item When the variable is decided in $t-1$, such as the capital stock in our simple model, we write ${\tt x(-1)}$. \item Finally, when a variable is decided in the next period, $t+1$, such as consumption in the Euler equation, we write ${\tt x(+1)}$.
 \end{itemize}
 Hence the required code to declare our model in {\sc Dynare} will be:
@@ -164,7 +164,7 @@ end;
 \end{verbatim}
 so that the level of consumption is actually given by ${\tt exp(c)}$.
 \paragraph{Solving the model}
-\begin{enumerate} 
+\begin{enumerate}
 \item Now we need to provide numerical initial conditions for the computation of the deterministic steady state. This is done with the sequence between {\tt initval;} and {\tt end;}. Each variable, endogenous or exogenous, should be initialized. In our example, we give the exact values of the deterministic equilibrium in absence of shocks. This takes the form
 \begin{verbatim}
 initval;
@@ -206,12 +206,12 @@ Number of periods on which to compute the IRFs (default = 40)
 \item {\tt nofunctions}:
 Doesn't print the coefficients of the approximated solution
 \item {\tt nomoments}:
-Doesn't print moments of the endogenous variables 
+Doesn't print moments of the endogenous variables
 \item {\tt order} = [1,2,3]:
 Order of Taylor approximation (default = 2)
 \item {\tt replic} = Integer:
 Number of simulated series used to compute the IRFs (default = 1, if order = 1, and 50 otherwise)
-\end{itemize} 
+\end{itemize}
 
 In our first example, we use simply:
 \begin{verbatim}
@@ -324,8 +324,8 @@ end;
 
 stoch_simul(periods=2000, drop=200);
 \end{verbatim}
-\bibliographystyle{Usmad}
-\bibliography{/papers/biblio/michel}
+%\bibliographystyle{Usmad}
+%\bibliography{/papers/biblio/michel}
 
 \end{document}
 

doc/manual/source/bibliography.rst

@@ -8,10 +8,12 @@ Bibliography
 * Adjemian, Stéphane, Matthieu Darracq Parriès and Stéphane Moyen (2008): “Towards a monetary policy evaluation framework”, *European Central Bank Working Paper*, 942.
 * Aguiar, Mark and Gopinath, Gita (2004): “Emerging Market Business Cycles: The Cycle is the Trend,” *NBER* Working Paper, 10734.
 * Amisano, Gianni and Tristani, Oreste (2010): “Euro area inflation persistence in an estimated nonlinear DSGE model”, *Journal of Economic Dynamics and Control*, 34(10), 1837–1858.
-* Andreasen, Martin M., Jesús Fernández-Villaverde, and Juan Rubio-Ramírez (2018): “The Pruned State-Space System for Non-Linear DSGE Models: Theory and Empirical Applications,” *Review of Economic Studies*, 85(1), pp. 1-49.
+* Andreasen, Martin M., Jesús Fernández-Villaverde, and Juan Rubio-Ramírez (2018): “The Pruned State-Space System for Non-Linear DSGE Models: Theory and Empirical Applications,” *Review of Economic Studies*, 85(1), 1-49.
+* Andrle, Michal and Miroslav Plašil (2018): “Econometrics with system priors,” *Economics Letters*, 172, 134-137.
 * Andrews, Donald W.K (1991): “Heteroskedasticity and autocorrelation consistent covariance matrix estimation”, *Econometrica*, 59(3), 817–858.
 * Backus, David K., Patrick J. Kehoe, and Finn E. Kydland (1992): “International Real Business Cycles,” *Journal of Political Economy*, 100(4), 745–775.
 * Baxter, Marianne and Robert G. King (1999): “Measuring Business Cycles: Approximate Band-pass Filters for Economic Time Series,” *Review of Economics and Statistics*, 81(4), 575–593.
+* Born, Benjamin and Johannes Pfeifer (2014): “Policy risk and the business cycle”, *Journal of Monetary Economics*, 68, 68-85.
 * Boucekkine, Raouf (1995): “An alternative methodology for solving nonlinear forward-looking models,” *Journal of Economic Dynamics and Control*, 19, 711–734.
 * Brooks, Stephen P., and Andrew Gelman (1998): “General methods for monitoring convergence of iterative simulations,” *Journal of Computational and Graphical Statistics*, 7, pp. 434–455.
 * Cardoso, Margarida F., R. L. Salcedo and S. Feyo de Azevedo (1996): “The simplex simulated annealing approach to continuous non-linear optimization,” *Computers & Chemical Engineering*, 20(9), 1065-1080.
@@ -21,19 +23,25 @@ Bibliography
 * Collard, Fabrice (2001): “Stochastic simulations with Dynare: A practical guide”.
 * Collard, Fabrice and Michel Juillard (2001a): “Accuracy of stochastic perturbation methods: The case of asset pricing models,” *Journal of Economic Dynamics and Control*, 25, 979–999.
 * Collard, Fabrice and Michel Juillard (2001b): “A Higher-Order Taylor Expansion Approach to Simulation of Stochastic Forward-Looking Models with an Application to a Non-Linear Phillips Curve,” *Computational Economics*, 17, 125–139.
-* Corona, Angelo, M. Marchesi, Claudio Martini, and Sandro Ridella (1987): “Minimizing multimodal functions of continuous variables with the “simulated annealing” algorithm”, *ACM Transactions on Mathematical Software*, 13(3), 262–280.
-* Del Negro, Marco and Franck Schorfheide (2004): “Priors from General Equilibrium Models for VARs”, *International Economic Review*, 45(2), 643–673.
+* Corana, Angelo, M. Marchesi, Claudio Martini, and Sandro Ridella (1987): “Minimizing multimodal functions of continuous variables with the “simulated annealing” algorithm”, *ACM Transactions on Mathematical Software*, 13(3), 262–280.
+* Cuba-Borda, Pablo, Luca Guerrieri, and Matteo Iacoviello (2019): "Likelihood evaluation of models with occasionally binding constraints", Journal of Applied Econometrics, 34(7), 1073-1085
+* Del Negro, Marco and Frank Schorfheide (2004): “Priors from General Equilibrium Models for VARs”, *International Economic Review*, 45(2), 643–673.
 * Dennis, Richard (2007): “Optimal Policy In Rational Expectations Models: New Solution Algorithms”, *Macroeconomic Dynamics*, 11(1), 31–55.
+* Duffie, Darrel and Kenneth J. Singleton (1993): “Simulated Moments Estimation of Markov Models of Asset Prices”, *Econometrica*, 61(4), 929-952.
 * Durbin, J. and S. J. Koopman (2012), *Time Series Analysis by State Space Methods*, Second Revised Edition, Oxford University Press.
 * Fair, Ray and John Taylor (1983): “Solution and Maximum Likelihood Estimation of Dynamic Nonlinear Rational Expectation Models,” *Econometrica*, 51, 1169–1185.
+* Fernández-Villaverde, Jesús (2010): “The econometrics of DSGE models,” *SERIEs*, 1, 3–49.
 * Fernández-Villaverde, Jesús and Juan Rubio-Ramírez (2004): “Comparing Dynamic Equilibrium Economies to Data: A Bayesian Approach,” *Journal of Econometrics*, 123, 153–187.
 * Fernández-Villaverde, Jesús and Juan Rubio-Ramírez (2005): “Estimating Dynamic Equilibrium Economies: Linear versus Nonlinear Likelihood,” *Journal of Applied Econometrics*, 20, 891–910.
-* Fernández-Villaverde, Jesús (2010): “The econometrics of DSGE models,” *SERIEs*, 1, 3–49.
 * Ferris, Michael C. and Todd S. Munson (1999): “Interfaces to PATH 3.0: Design, Implementation and Usage”, *Computational Optimization and Applications*, 12(1), 207–227.
+* Galí, Jordi (2015): “Monetary Policy, Inflation, and the Business Cycle,” 2nd Edition, Princeton University Press, Princeton.
 * Geweke, John (1992): “Evaluating the accuracy of sampling-based approaches to the calculation of posterior moments,” in J.O. Berger, J.M. Bernardo, A.P. Dawid, and A.F.M. Smith (eds.) Proceedings of the Fourth Valencia International Meeting on Bayesian Statistics, pp. 169–194, Oxford University Press.
 * Geweke, John (1999): “Using simulation methods for Bayesian econometric models: Inference, development and communication,” *Econometric Reviews*, 18(1), 1–73.
+* Giovannini, Massimo, Philipp Pfeiffer, and Marco Ratto (2021), “Efficient and robust inference of models with occasionally binding constraints,” Working Papers 2021-03, Joint Research Centre, European Commission
 * Giordani, Paolo, Michael Pitt, and Robert Kohn (2011): “Bayesian Inference for Time Series State Space Models” in: *The Oxford Handbook of Bayesian Econometrics*, ed. by John Geweke, Gary Koop, and Herman van Dijk, Oxford University Press, 61–124.
+* Guerrieri, Luca and Matteo Iacoviello (2015): “OccBin: A toolkit for solving dynamic models with occasionally binding constraints easily,” *Journal of Monetary Economics*, 70, 22–38.
 * Goffe, William L., Gary D. Ferrier, and John Rogers (1994): “Global Optimization of Statistical Functions with Simulated Annealing,” *Journal of Econometrics*, 60(1/2), 65–100.
+* Hansen, Lars P. (1982): “Large sample properties of generalized method of moments estimators,” Econometrica, 50(4), 1029–1054.
 * Hansen, Nikolaus and Stefan Kern (2004): “Evaluating the CMA Evolution Strategy on Multimodal Test Functions”. In: *Eighth International Conference on Parallel Problem Solving from Nature PPSN VIII*, Proceedings, Berlin: Springer, 282–291.
 * Harvey, Andrew C. and Garry D.A. Phillips (1979): “Maximum likelihood estimation of regression models with autoregressive-moving average disturbances,” *Biometrika*, 66(1), 49–58.
 * Herbst, Edward (2015): “Using the “Chandrasekhar Recursions” for Likelihood Evaluation of DSGE Models,” *Computational Economics*, 45(4), 693–705.
@@ -41,8 +49,8 @@ Bibliography
 * Iskrev, Nikolay (2010): “Local identification in DSGE models,” *Journal of Monetary Economics*, 57(2), 189–202.
 * Judd, Kenneth (1996): “Approximation, Perturbation, and Projection Methods in Economic Analysis”, in *Handbook of Computational Economics*, ed. by Hans Amman, David Kendrick, and John Rust, North Holland Press, 511–585.
 * Juillard, Michel (1996): “Dynare: A program for the resolution and simulation of dynamic models with forward variables through the use of a relaxation algorithm,” CEPREMAP, *Couverture Orange*, 9602.
-* Kim, Jinill and Sunghyun Kim (2003): “Spurious welfare reversals in international business cycle models,” *Journal of International Economics*, 60, 471–500.
 * Kanzow, Christian and Stefania Petra (2004): “On a semismooth least squares formulation of complementarity problems with gap reduction,” *Optimization Methods and Software*, 19, 507–525.
+* Kim, Jinill and Sunghyun Kim (2003): “Spurious welfare reversals in international business cycle models,” *Journal of International Economics*, 60, 471–500.
 * Kim, Jinill, Sunghyun Kim, Ernst Schaumburg, and Christopher A. Sims (2008): “Calculating and using second-order accurate solutions of discrete time dynamic equilibrium models,” *Journal of Economic Dynamics and Control*, 32(11), 3397–3414.
 * Komunjer, Ivana and Ng, Serena (2011): ”Dynamic identification of dynamic stochastic general equilibrium models”, *Econometrica*, 79, 1995–2032.
 * Koop, Gary (2003), *Bayesian Econometrics*, John Wiley & Sons.
@@ -51,25 +59,26 @@ Bibliography
 * Kuntsevich, Alexei V. and Franz Kappel (1997): “SolvOpt - The solver for local nonlinear optimization problems (version 1.1, Matlab, C, FORTRAN)”, University of Graz, Graz, Austria.
 * Laffargue, Jean-Pierre (1990): “Résolution d’un modèle macroéconomique avec anticipations rationnelles”, *Annales d’Économie et Statistique*, 17, 97–119.
 * Liu, Jane and Mike West (2001): “Combined parameter and state estimation in simulation-based filtering”, in *Sequential Monte Carlo Methods in Practice*, Eds. Doucet, Freitas and Gordon, Springer Verlag.
-* Lubik, Thomas and Frank Schorfheide (2007): “Do Central Banks Respond to Exchange Rate Movements? A Structural Investigation,” *Journal of Monetary Economics*, 54(4), 1069–1087.
 * Murray, Lawrence M., Emlyn M. Jones and John Parslow (2013): “On Disturbance State-Space Models and the Particle Marginal Metropolis-Hastings Sampler”, *SIAM/ASA Journal on Uncertainty Quantification*, 1, 494–521.
 * Mutschler, Willi (2015): “Identification of DSGE models - The effect of higher-order approximation and pruning“, *Journal of Economic Dynamics & Control*, 56, 34-54.
+* Mutschler, Willi (2018): “Higher-order statistics for DSGE models”, *Econometrics and Statistics*, 6(C), 44-56.
 * Pearlman, Joseph, David Currie, and Paul Levine (1986): “Rational expectations models with partial information,” *Economic Modelling*, 3(2), 90–105.
 * Planas, Christophe, Marco Ratto and Alessandro Rossi (2015): “Slice sampling in Bayesian estimation of DSGE models”.
 * Pfeifer, Johannes (2013): “A Guide to Specifying Observation Equations for the Estimation of DSGE Models”.
 * Pfeifer, Johannes (2014): “An Introduction to Graphs in Dynare”.
 * Qu, Zhongjun and Tkachenko, Denis (2012): “Identification and frequency domain quasi-maximum likelihood estimation of linearized dynamic stochastic general equilibrium models“, *Quantitative Economics*, 3, 95–132.
-* Rabanal, Pau and Juan Rubio-Ramirez (2003): “Comparing New Keynesian Models of the Business Cycle: A Bayesian Approach,” Federal Reserve of Atlanta, *Working Paper Series*, 2003-30.
+* Rabanal, Pau and Juan Rubio-Ramírez (2003): “Comparing New Keynesian Models of the Business Cycle: A Bayesian Approach,” Federal Reserve of Atlanta, *Working Paper Series*, 2003-30.
 * Raftery, Adrian E. and Steven Lewis (1992): “How many iterations in the Gibbs sampler?,” in *Bayesian Statistics, Vol. 4*, ed. J.O. Berger, J.M. Bernardo, A.P. * Dawid, and A.F.M. Smith, Clarendon Press: Oxford, pp. 763-773.
 * Ratto, Marco (2008): “Analysing DSGE models with global sensitivity analysis”, *Computational Economics*, 31, 115–139.
 * Ratto, Marco and Iskrev, Nikolay (2011): “Identification Analysis of DSGE Models with DYNARE.“, *MONFISPOL* 225149.
-* Schorfheide, Frank (2000): “Loss Function-based evaluation of DSGE models,” *Journal of Applied Econometrics*, 15(6), 645–670.
+* Ruge-Murcia, Francisco J. (2012): “Estimating nonlinear DSGE models by the simulated method of moments: With an application to business cycles“, *Journal of Economic Dynamics and Control*, 36, 914-938.
 * Schmitt-Grohé, Stephanie and Martin Uríbe (2004): “Solving Dynamic General Equilibrium Models Using a Second-Order Approximation to the Policy Function,” *Journal of Economic Dynamics and Control*, 28(4), 755–775.
 * Schnabel, Robert B. and Elizabeth Eskow (1990): “A new modified Cholesky algorithm,” *SIAM Journal of Scientific and Statistical Computing*, 11, 1136–1158.
+* Schorfheide, Frank (2000): “Loss Function-based evaluation of DSGE models,” *Journal of Applied Econometrics*, 15(6), 645–670.
 * Sims, Christopher A., Daniel F. Waggoner and Tao Zha (2008): “Methods for inference in large multiple-equation Markov-switching models,” *Journal of Econometrics*, 146, 255–274.
 * Skoeld, Martin and Gareth O. Roberts (2003): “Density Estimation for the Metropolis-Hastings Algorithm,” *Scandinavian Journal of Statistics*, 30, 699–718.
 * Smets, Frank and Rafael Wouters (2003): “An Estimated Dynamic Stochastic General Equilibrium Model of the Euro Area,” *Journal of the European Economic Association*, 1(5), 1123–1175.
 * Stock, James H. and Mark W. Watson (1999). “Forecasting Inflation,”, *Journal of Monetary Economics*, 44(2), 293–335.
 * Uhlig, Harald (2001): “A Toolkit for Analysing Nonlinear Dynamic Stochastic Models Easily,” in *Computational Methods for the Study of Dynamic Economies*, Eds. Ramon Marimon and Andrew Scott, Oxford University Press, 30–61.
-* U.S. Census Bureau (2017): “X-13 ARIMA-SEATSReference Manual”.
+* U.S. Census Bureau (2017): “X-13 ARIMA-SEATS Reference Manual”.
 * Villemot, Sébastien (2011): “Solving rational expectations models at first order: what Dynare does,” *Dynare Working Papers*, 2, CEPREMAP.

doc/manual/source/examples.rst

@@ -51,12 +51,23 @@ description, please refer to the comments inside the files themselves.
     Small open economy RBC model with shocks to the growth trend,
     presented in *Aguiar and Gopinath (2004)*.
 
+``Gali_2015.mod``
+
+    Basic New Keynesian model of *Galí (2015)*, Chapter 3 showing how to
+    i) use "system prior"-type prior restrictions as in *Andrle and Plašil (2018)*
+    and ii) run prior/posterior-functions.
+
 ``NK_baseline.mod``
 
     Baseline New Keynesian Model estimated in *Fernández-Villaverde
     (2010)*. It demonstrates how to use an explicit steady state file
     to update parameters and call a numerical solver.
 
+``Occbin_example.mod``
+
+    RBC model with two occasionally binding constraints. Demonstrates
+    how to set up Occbin.
+
 ``Ramsey_Example.mod``
 
     File demonstrating how to conduct optimal policy experiments in a 
@@ -67,4 +78,4 @@ description, please refer to the comments inside the files themselves.
 
     File demonstrating how to conduct optimal policy experiments in a 
     simple New Keynesian model under commitment (Ramsey) with a user-defined
-    conditional steady state file
+    conditional steady state file

doc/manual/source/time-series.rst

@@ -15,7 +15,7 @@ class and methods for dates. Below, you will first find the class and
 methods used for creating and dealing with dates and then the class
 used for using time series. Dynare also provides an interface to the
 X-13 ARIMA-SEATS seasonal adjustment program produced, distributed, and
-maintained by the US Census Bureau.
+maintained by the US Census Bureau (2017).
 
 
 Dates

doc/manual/utils/dynare_lex.py

@@ -36,8 +36,8 @@ class DynareLexer(RegexLexer):
         "dynare","var","varexo","varexo_det","parameters","change_type","model_local_variable",
         "predetermined_variables","trend_var","log_trend_var","external_function",
         "write_latex_original_model","write_latex_dynamic_model",
-        "write_latex_static_model","resid","initval_file","histval_file","dsample",
-        "periods","values","corr","steady","check","model_diagnostics","model_info",
+        "write_latex_static_model","write_latex_steady_state_model","resid","initval_file","histval_file","dsample",
+        "periods","values","scales","corr","stderr","steady","check","model_diagnostics","model_info",
         "print_bytecode_dynamic_model"," print_bytecode_static_model",
         "perfect_foresight_setup","perfect_foresight_solver","simul","stoch_simul",
         "extended_path","varobs","estimation","unit_root_vars","bvar_density",
@@ -52,7 +52,8 @@ class DynareLexer(RegexLexer):
         "save_params_and_steady_state","load_params_and_steady_state",
         "dynare_version","write_latex_definitions","write_latex_parameter_table",
         "write_latex_prior_table","collect_latex_files","prior_function",
-        "posterior_function","generate_trace_plots","evaluate_planner_objective")
+        "posterior_function","generate_trace_plots","evaluate_planner_objective",
+        "occbin_setup","occbin_solver","occbin_write_regimes","occbin_graph","method_of_moments")
 
     report_commands = ("report","addPage","addSection","addGraph","addTable",
         "addSeries","addParagraph","addVspace","write","compile")
@@ -80,6 +81,7 @@ class DynareLexer(RegexLexer):
                 'shock_groups','conditional_forecast_paths','optim_weights',
                 'osr_params_bounds','ramsey_constraints','irf_calibration',
                 'moment_calibration','identification','svar_identification',
+                'matched_moments','occbin_constraints','surprise','overwrite','bind','relax',
                 'verbatim','end','node','cluster','paths','hooks'), prefix=r'\b', suffix=r'\s*\b'),Keyword.Reserved),
 
             # FIXME: Commands following multiline comments are not highlighted properly.

doc/parallel/parallel.tex

@@ -1,14 +1,15 @@
 % ----------------------------------------------------------------
 % AMS-LaTeX Paper ************************************************
 % **** -----------------------------------------------------------
-\documentclass[12pt,a4paper,pdftex,nofootinbib]{article}
-\usepackage[cp1252]{inputenc}
+\documentclass[12pt,a4paper,pdftex]{article}
+\usepackage[margin=2.5cm]{geometry}
+\usepackage[utf8]{inputenc}
 \usepackage{amssymb,amsmath}
-\usepackage[pdftex]{graphicx}
+\usepackage{graphicx}
 \usepackage{epstopdf}
 \usepackage{natbib}
 \usepackage{verbatim}
-\usepackage[pdftex]{color}
+\usepackage{xcolor}
 \usepackage{psfrag}
 \usepackage{setspace}
 \usepackage{rotating}
@@ -49,6 +50,15 @@
 \def \supp{{\rm supp}}
 \def \var{{\rm var}}
 
+\usepackage[pdfpagelabels]{hyperref}
+\hypersetup{
+pdfproducer = {LaTeX},
+colorlinks,
+linkcolor=blue,
+filecolor=yellow,
+urlcolor=green,
+citecolor=green}
+
 % ----------------------------------------------------------------
 \begin{document}
 
@@ -349,7 +359,7 @@ Finally, the DYNARE command line options are:
   \item \verb"parallel": trigger the parallel computation using the first cluster specified in config file
   \item \verb"parallel=<clustername>": trigger the parallel computation, using the given cluster
   \item \verb"parallel_slave_open_mode": use the leaveSlaveOpen mode in the cluster
-  \item \verb"parallel_test": just test the cluster, don’t actually run the MOD file
+  \item \verb"parallel_test": just test the cluster, don't actually run the MOD file
 
  \end{itemize}
 
@@ -828,7 +838,7 @@ The modified \verb"random_walk_metropolis_hastings.m" is therefore:
 \noindent\begin{tabular}[b]{| p{\linewidth} |}
 \hline
 \begin{verbatim}
-function random_walk_metropolis_hastings(TargetFun,ProposalFun,…,varargin)
+function random_walk_metropolis_hastings(TargetFun,ProposalFun,\ldots,varargin)
 [...]
 % here we wrap all local variables needed by the <*>_core function
 localVars = struct('TargetFun', TargetFun, ...
@@ -970,11 +980,11 @@ On the other hand, under the parallel implementation, a parallel monitoring plot
 
 \section{Parallel DYNARE: testing}
 We checked the new parallel platform for DYNARE performing a number of tests, using different models and computer architectures. We present here all tests performed with Windows XP/MATLAB. However, similar tests were performed successfully under Linux/Ubuntu environment.
-In the Bayesian estimation of DSGE models with DYNARE, most of the computing time is devoted to the posterior parameter estimation with the Metropolis algorithm. The first and second tests are therefore focused on the parallelization of the Random Walking Metropolis Hastings algorithm (Sections \ref{s:test1}-\ref{s:test2}). In addition, further tests (Sections \ref{s:test3}-\ref{s:test4}) are devoted to test all the parallelized functions in DYNARE. Finally, we compare the two parallel implementations of the Metropolis Hastings algorithms, available in DYNARE: the Independent and the Random Walk (Section \ref{s:test5}).
+In the Bayesian estimation of DSGE models with DYNARE, most of the computing time is devoted to the posterior parameter estimation with the Metropolis algorithm. The first and second tests are therefore focused on the parallelization of the Random Walking Metropolis Hastings algorithm (Sections \ref{s:test1}-\ref{s:test2}). In addition, further tests (Sections \ref{s:test3}-\ref{s:test4}) are devoted to test all the parallelized functions in DYNARE. %Finally, we compare the two parallel implementations of the Metropolis Hastings algorithms, available in DYNARE: the Independent and the Random Walk (Section \ref{s:test5}).
 
 \subsection{Test 1.}\label{s:test1}
 The main goal here was to evaluate the parallel package on a \emph{fixed hardware platform} and using chains of \emph{variable length}. The model used for testing is a modification of \cite{Hradisky_etal_2006}. This is a small scale open economy DSGE model with 6 observed variables, 6 endogenous variables and 19 parameters to be estimated.
-We estimated the model on a bi-processor machine (Fujitsu Siemens, Celsius R630) powered with an Intel(R) Xeon(TM) CPU 2.80GHz Hyper Treading Technology; first with the original serial Metropolis and subsequently using the parallel solution, to take advantage of the two processors technology. We ran chains of increasing length: 2500, 5000, 10,000, 50,000, 100,000, 250,000, 1,000,000.
+We estimated the model on a bi-processor machine (Fujitsu Siemens, Celsius R630) powered with an Intel\textsuperscript{\textregistered} Xeon\texttrademark CPU 2.80GHz Hyper Treading Technology; first with the original serial Metropolis and subsequently using the parallel solution, to take advantage of the two processors technology. We ran chains of increasing length: 2500, 5000, 10,000, 50,000, 100,000, 250,000, 1,000,000.
 
 \begin{figure}[!ht]
 \begin{centering}
@@ -997,8 +1007,8 @@ Overall results are given in Figure \ref{fig:test_time_comp}, showing the comput
 The scope of the second test was to verify if results were robust over different hardware platforms.
 We estimated the model with chain lengths of 1,000,000 runs on the following hardware platforms:
 \begin{itemize}
-\item Single processor machine: Intel(R) Pentium4(R) CPU 3.40GHz with Hyper Treading Technology (Fujitsu-Siemens Scenic Esprimo);
-\item Bi-processor machine: two CPU's Intel(R) Xeon(TM) 2.80GHz Hyper Treading Technology (Fujitsu-Siemens, Celsius R630);
+\item Single processor machine: Intel\textsuperscript{\textregistered} Pentium4\textsuperscript{\textregistered} CPU 3.40GHz with Hyper Treading Technology (Fujitsu-Siemens Scenic Esprimo);
+\item Bi-processor machine: two CPU's Intel\textsuperscript{\textregistered} Xeon\texttrademark 2.80GHz Hyper Treading Technology (Fujitsu-Siemens, Celsius R630);
 \item Dual core machine: Intel Centrino T2500 2.00GHz Dual Core  (Fujitsu-Siemens, LifeBook S Series).
 \end{itemize}
 
@@ -1042,7 +1052,7 @@ Unplugged network cable. &
 Given the excellent results reported above, we have parallelized many other DYNARE functions. This implies that parallel instances can be invoked many times during a single DYNARE session. Under the basic parallel toolbox implementation, that we call the `Open/Close' strategy, this implies that MATLAB instances are opened and closed many times by system calls, possibly slowing down the computation, specially for `entry-level' computer resources. As mentioned before, this suggested to implement an alternative strategy for the parallel toolbox, that we call the `Always-Open' strategy, where the slave MATLAB threads, once opened, stay alive and wait for new tasks assigned by the master until the full DYNARE procedure is completed. We show next the tests of these latest implementations.
 
 \subsection{Test 3}\label{s:test3}
-In this Section we use the \cite{Lubik2003} model as test function\footnote{The \cite{Lubik2003} model is also selected as the `official' test model for the parallel toolbox in DYNARE.} and a very simple computer class, quite diffuse nowadays: Netbook personal Computer. In particular we used the Dell Mini 10 with Processor Intel® Atom™ Z520 (1,33 GHz, 533 MHz), 1 GB di RAM (with Hyper-trading). First, we tested the computational gain of running a full Bayesian estimation: Metropolis (two parallel chains), MCMC diagnostics, posterior IRF's and filtered, smoothed, forecasts, etc. In other words, we designed DYNARE sessions that invoke all parallelized functions. Results are shown in Figures \ref{fig:netbook_complete_openclose}-\ref{fig:netbook_partial_openclose}.
+In this Section we use the \cite{Lubik2003} model as test function\footnote{The \cite{Lubik2003} model is also selected as the `official' test model for the parallel toolbox in DYNARE.} and a very simple computer class, quite diffuse nowadays: Netbook personal Computer. In particular we used the Dell Mini 10 with Processor Intel\textsuperscript{\textregistered} Atom\texttrademark Z520 (1,33 GHz, 533 MHz), 1 GB di RAM (with Hyper-trading). First, we tested the computational gain of running a full Bayesian estimation: Metropolis (two parallel chains), MCMC diagnostics, posterior IRF's and filtered, smoothed, forecasts, etc. In other words, we designed DYNARE sessions that invoke all parallelized functions. Results are shown in Figures \ref{fig:netbook_complete_openclose}-\ref{fig:netbook_partial_openclose}.
 \begin{figure}[p]
 \begin{centering}
   % Requires \usepackage{graphicx}
@@ -1143,49 +1153,49 @@ The methodology identified for parallelizing MATLAB codes within DYNARE proved t
 \begin{figure}
 \begin{centering}
   % Requires \usepackage{graphicx}
-  \epsfxsize=250pt \epsfbox{RWMH_quest1_PriorsAndPosteriors1Comp.pdf}
+  \epsfxsize=300pt \epsfbox{RWMH_quest1_PriorsAndPosteriors1Comp.pdf}
   \caption{Prior (grey lines) and posterior density of estimated parameters (black = 100,000 runs; red = 1,000,000 runs) using the RWMH algorithm \citep[QUEST III model][]{Ratto_et_al_EconModel2009}.}\label{fig:quest_RWMH_comp1}
 \end{centering}
 \end{figure}
 \begin{figure}
 \begin{centering}
   % Requires \usepackage{graphicx}
-  \epsfxsize=250pt \epsfbox{RWMH_quest1_PriorsAndPosteriors2Comp.pdf}
+  \epsfxsize=300pt \epsfbox{RWMH_quest1_PriorsAndPosteriors2Comp.pdf}
   \caption{Prior (grey lines) and posterior density of estimated parameters (black = 100,000 runs; red = 1,000,000 runs) using the RWMH algorithm \citep[QUEST III model][]{Ratto_et_al_EconModel2009}.}\label{fig:quest_RWMH_comp2}
 \end{centering}
 \end{figure}
 \begin{figure}
 \begin{centering}
   % Requires \usepackage{graphicx}
-  \epsfxsize=250pt \epsfbox{RWMH_quest1_PriorsAndPosteriors3Comp.pdf}
+  \epsfxsize=300pt \epsfbox{RWMH_quest1_PriorsAndPosteriors3Comp.pdf}
   \caption{Prior (grey lines) and posterior density of estimated parameters (black = 100,000 runs; red = 1,000,000 runs) using the RWMH algorithm \citep[QUEST III model][]{Ratto_et_al_EconModel2009}.}\label{fig:quest_RWMH_comp3}
 \end{centering}
 \end{figure}
 \begin{figure}
 \begin{centering}
   % Requires \usepackage{graphicx}
-  \epsfxsize=250pt \epsfbox{RWMH_quest1_PriorsAndPosteriors4Comp.pdf}
+  \epsfxsize=300pt \epsfbox{RWMH_quest1_PriorsAndPosteriors4Comp.pdf}
   \caption{Prior (grey lines) and posterior density of estimated parameters (black = 100,000 runs; red = 1,000,000 runs) using the RWMH algorithm \citep[QUEST III model][]{Ratto_et_al_EconModel2009}.}\label{fig:quest_RWMH_comp4}
 \end{centering}
 \end{figure}
 \begin{figure}
 \begin{centering}
   % Requires \usepackage{graphicx}
-  \epsfxsize=250pt \epsfbox{RWMH_quest1_PriorsAndPosteriors5Comp.pdf}
+  \epsfxsize=300pt \epsfbox{RWMH_quest1_PriorsAndPosteriors5Comp.pdf}
   \caption{Prior (grey lines) and posterior density of estimated parameters (black = 100,000 runs; red = 1,000,000 runs) using the RWMH algorithm \citep[QUEST III model][]{Ratto_et_al_EconModel2009}.}\label{fig:quest_RWMH_comp5}
 \end{centering}
 \end{figure}
 \begin{figure}
 \begin{centering}
   % Requires \usepackage{graphicx}
-  \epsfxsize=250pt \epsfbox{RWMH_quest1_PriorsAndPosteriors6Comp.pdf}
+  \epsfxsize=300pt \epsfbox{RWMH_quest1_PriorsAndPosteriors6Comp.pdf}
   \caption{Prior (grey lines) and posterior density of estimated parameters (black = 100,000 runs; red = 1,000,000 runs) using the RWMH algorithm \citep[QUEST III model][]{Ratto_et_al_EconModel2009}.}\label{fig:quest_RWMH_comp6}
 \end{centering}
 \end{figure}
 \begin{figure}
 \begin{centering}
   % Requires \usepackage{graphicx}
-  \epsfxsize=250pt \epsfbox{RWMH_quest1_PriorsAndPosteriors7Comp.pdf}
+  \epsfxsize=300pt \epsfbox{RWMH_quest1_PriorsAndPosteriors7Comp.pdf}
   \caption{Prior (grey lines) and posterior density of estimated parameters (black = 100,000 runs; red = 1,000,000 runs) using the RWMH algorithm \citep[QUEST III model][]{Ratto_et_al_EconModel2009}.}\label{fig:quest_RWMH_comp7}
 \end{centering}
 \end{figure}

examples/Gali_2015.mod

@@ -0,0 +1,184 @@
+/*
+ * This file shows how to use "system prior"-type prior restrictions as in
+ * Michal Andrle/Miroslav Plašil (2018): "Econometrics with system priors",
+ * Economics Letters, 172, pp. 134-137 during estimation based on
+ * the baseline New Keynesian model of Jordi Galí (2015): Monetary Policy, Inflation,
+ * and the Business Cycle, Princeton University Press, Second Edition, Chapter 3
+ *
+ * THIS MOD-FILE REQUIRES DYNARE 4.5 OR HIGHER
+ *
+ * Notes:
+ *  - The estimation will automatically take the Gali_2015_prior_restrictions.m into
+ *      account, which has the required name and format
+ *  - Estimation is based on simulated data
+ *  - The file also shows how to use a prior/posterior-function
+ *
+ * This implementation was written by Johannes Pfeifer. In case you spot mistakes,
+ * email me at jpfeifer@gmx.de
+ *
+ * Please note that the following copyright notice only applies to this Dynare
+ */
+
+/*
+ * Copyright (C) 2021 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/>.
+ */
+
+var pi          ${\pi}$                 (long_name='inflation')
+    y_gap       ${\tilde y}$            (long_name='output gap')
+    y_nat       ${y^{nat}}$             (long_name='natural output')      //(in contrast to the textbook defined in deviation from steady state)
+    y           ${y}$                   (long_name='output')
+    yhat        ${\hat y}$              (long_name='output deviation from steady state')
+    r_nat       ${r^{nat}}$             (long_name='natural interest rate')
+    i           ${i}$                   (long_name='nominal interrst rate')
+    n           ${n}$                   (long_name='hours worked')
+    nu          ${\nu}$                 (long_name='AR(1) monetary policy shock process')
+    a           ${a}$                   (long_name='AR(1) technology shock process')
+    z           ${z}$                   (long_name='AR(1) preference shock process')
+    p           ${p}$                   (long_name='price level')
+    w           ${w}$                   (long_name='nominal wage')
+    c           ${c}$                   (long_name='consumption')
+;
+
+varexo  eps_a       ${\varepsilon_a}$       (long_name='technology shock')
+        eps_nu  ${\varepsilon_\nu}$     (long_name='monetary policy shock')
+        eps_z       ${\varepsilon_z}$   (long_name='preference shock innovation')
+       ;
+
+parameters alppha       ${\alpha}$     (long_name='capital share')
+    betta               ${\beta}$       (long_name='discount factor')
+    rho_a               ${\rho_a}$      (long_name='autocorrelation technology shock')
+    rho_nu          ${\rho_{\nu}}$  (long_name='autocorrelation monetary policy shock')
+    rho_z               ${\rho_{z}}$    (long_name='autocorrelation monetary demand shock')
+    siggma              ${\sigma}$      (long_name='inverse EIS')
+    varphi              ${\varphi}$     (long_name='inverse Frisch elasticity')
+    phi_pi              ${\phi_{\pi}}$  (long_name='inflation feedback Taylor Rule')
+    phi_y               ${\phi_{y}}$    (long_name='output feedback Taylor Rule')
+    eta                 ${\eta}$        (long_name='semi-elasticity of money demand')
+    epsilon             ${\epsilon}$    (long_name='demand elasticity')
+    theta               ${\theta}$      (long_name='Calvo parameter')
+    ;
+%----------------------------------------------------------------
+% Parametrization, p. 67  and p. 113-115
+%----------------------------------------------------------------
+siggma = 1;
+varphi = 5;
+phi_pi = 1.5;
+phi_y  = 0.125;
+theta = 3/4;
+rho_nu =0.5;
+rho_z  = 0.5;
+rho_a  = 0.9;
+betta  = 0.99;
+eta  = 3.77; %footnote 11, p. 115
+alppha = 1/4;
+epsilon = 9;
+
+%----------------------------------------------------------------
+% First Order Conditions
+%----------------------------------------------------------------
+
+model(linear);
+//Composite parameters
+#Omega=(1-alppha)/(1-alppha+alppha*epsilon);        %defined on page 60
+#psi_n_ya=(1+varphi)/(siggma*(1-alppha)+varphi+alppha);   %defined on page 62
+#lambda=(1-theta)*(1-betta*theta)/theta*Omega;      %defined on page 61
+#kappa=lambda*(siggma+(varphi+alppha)/(1-alppha));     %defined on page 63
+[name='New Keynesian Phillips Curve eq. (22)']
+pi=betta*pi(+1)+kappa*y_gap;
+[name='Dynamic IS Curve eq. (23)']
+y_gap=-1/siggma*(i-pi(+1)-r_nat)+y_gap(+1);
+[name='Interest Rate Rule eq. (26)']
+i=phi_pi*pi+phi_y*yhat+nu;
+[name='Definition natural rate of interest eq. (24)']
+r_nat=-siggma*psi_n_ya*(1-rho_a)*a+(1-rho_z)*z;
+[name='Definition natural output, eq. (20)']
+y_nat=psi_n_ya*a;
+[name='Definition output gap']
+y_gap=y-y_nat;
+[name='Monetary policy shock']
+nu=rho_nu*nu(-1)+eps_nu;
+[name='TFP shock']
+a=rho_a*a(-1)+eps_a;
+[name='Production function (eq. 14)']
+y=a+(1-alppha)*n;
+[name='Preference shock, p. 54']
+z=rho_z*z(-1) - eps_z;
+[name='Output deviation from steady state']
+yhat=y-steady_state(y);
+[name='Definition price level']
+pi=p-p(-1);
+[name='resource constraint, eq. (12)']
+y=c;
+[name='FOC labor, eq. (2)']
+w-p=siggma*c+varphi*n;
+end;
+
+
+shocks;
+    var eps_nu = 0.0025^2; //1 standard deviation shock of 25 basis points, i.e. 1 percentage point annualized
+end;
+
+% simulate data
+stoch_simul(periods=100,drop=0,irf=0) yhat;
+% save data
+datatomfile('sim_data',{'yhat'});
+
+estimated_params;
+theta,0.75,beta_pdf,0.5,0.1;
+betta, beta_pdf, 0.993, 0.002;
+alppha, beta_pdf, 0.25, 0.02;
+end;
+
+varobs yhat;
+
+% Run prior function to get prior slope of the PC based on independent priors
+hh=figure('Name','Slope of the Phillips Curve');
+prior_function(function='Gali_2015_PC_slope');
+PC_slope_vec=cell2mat(oo_.prior_function_results(:,1));
+
+optimal_bandwidth = mh_optimal_bandwidth(PC_slope_vec,length(PC_slope_vec),0,'gaussian');
+[density(:,1),density(:,2)] = kernel_density_estimate(PC_slope_vec,512,length(PC_slope_vec),optimal_bandwidth,'gaussian');
+figure(hh)
+subplot(3,1,1)
+plot(density(:,1),density(:,2));
+title('Prior')
+
+% Run estimation with 1 observation to show effect of _prior_restriction .m
+% on independent prior
+estimation(datafile='sim_data',mode_compute=5,mh_replic=2001,mh_nblocks=1,diffuse_filter,nobs=1,mh_jscale=0.8);
+posterior_function(function='Gali_2015_PC_slope');
+PC_slope_vec=cell2mat(oo_.posterior_function_results(:,1));
+optimal_bandwidth = mh_optimal_bandwidth(PC_slope_vec,length(PC_slope_vec),0,'gaussian');
+[density(:,1),density(:,2)] = kernel_density_estimate(PC_slope_vec,512,length(PC_slope_vec),optimal_bandwidth,'gaussian');
+figure(hh)
+subplot(3,1,2)
+plot(density(:,1),density(:,2));
+title('Updated Prior')
+
+
+% Run estimation with full observations
+estimation(datafile='sim_data',mode_compute=5,mh_replic=2001,mh_nblocks=1,diffuse_filter,nobs=100,mh_jscale=0.8);
+
+posterior_function(function='Gali_2015_PC_slope');
+PC_slope_vec=cell2mat(oo_.posterior_function_results(:,1));
+optimal_bandwidth = mh_optimal_bandwidth(PC_slope_vec,length(PC_slope_vec),0,'gaussian');
+[density(:,1),density(:,2)] = kernel_density_estimate(PC_slope_vec,512,length(PC_slope_vec),optimal_bandwidth,'gaussian');
+figure(hh)
+subplot(3,1,3)
+plot(density(:,1),density(:,2));
+title('Posterior')

examples/Gali_2015_PC_slope.m

@@ -0,0 +1,56 @@
+function output_cell =PC_slope(xparam1,M_,options_,oo_,estim_params_,bayestopt_,dataset_,dataset_info)
+% output_cell =PC_slope(xparam1,M_,options_,oo_,estim_params_,bayestopt_,dataset_,dataset_info);
+% This is an example file computing statistics on the prior/posterior draws. The
+% function allows read-only access to all Dynare structures. However, those
+% structures are local to this function.  Changing them will not affect
+% other Dynare functions and you cannot use them to pass results to other
+% Dynare functions.
+% The function takes one and only one output argument: an 1 by n cell.
+% Using functions like cell2mat, the contents of the cell can be easily
+% transformed back to matrices. See the fs2000_posterior_function.mod for
+% an example
+
+% INPUTS
+%   xparam1                      Current parameter draw
+%   M_           [structure]     Matlab's structure describing the Model (initialized by dynare, see @ref{M_}).
+%   options_     [structure]     Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
+%   oo_          [structure]     Matlab's structure gathering the results (initialized by dynare, see @ref{oo_}).
+%   estim_params_[structure]     Matlab's structure describing the estimated_parameters (initialized by dynare, see @ref{estim_params_}).
+%   bayestopt_   [structure]     Matlab's structure describing the parameter options (initialized by dynare, see @ref{bayestopt_}).
+%   dataset_     [structure]     Matlab's structure storing the dataset
+%   dataset_info [structure]     Matlab's structure storing the information about the dataset
+
+% Output
+%   output_cell  [1 by n cell]   1 by n Matlab cell allowing to store any
+%                                desired computation or result (strings, matrices, structures, etc.)
+
+% Copyright (C) 2021 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/>.
+
+
+%% store the slope based on the parameter draw
+NumberOfParameters = M_.param_nbr;
+for ii = 1:NumberOfParameters
+  paramname = deblank(M_.param_names{ii,:});
+  eval([ paramname ' = M_.params(' int2str(ii) ');']);
+end
+
+Omega=(1-alppha)/(1-alppha+alppha*epsilon);
+lambda=(1-theta)*(1-betta*theta)/theta*Omega;      %defined on page 61
+output_cell{1,1}=lambda*(siggma+(varphi+alppha)/(1-alppha));     %defined on page 63
+
+end

examples/Gali_2015_prior_restrictions.m

@@ -0,0 +1,40 @@
+function log_prior_val=Gali_2015_prior_restrictions(M_, oo_, options_, dataset_, dataset_info);
+% function prior_val=Gali_2015_prior_restrictions(M_, oo_, options_, dataset_, dataset_info);
+% Example of a _prior_restrictions-file automatically called during
+% estimation
+% It imposes a prior of the slope of the New Keynesian Phillips Curve of
+% 0.03. As the slope is a composite of other parameters with independent
+% priors, a separate function is required to do this.
+
+% Copyright (C) 2021 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 <https://www.gnu.org/licenses/>.
+
+
+% read out parameters to access them with their name
+NumberOfParameters = M_.param_nbr;
+for ii = 1:NumberOfParameters
+  paramname = M_.param_names{ii};
+  eval([ paramname ' = M_.params(' int2str(ii) ');']);
+end
+
+Omega=(1-alppha)/(1-alppha+alppha*epsilon);
+lambda=(1-theta)*(1-betta*theta)/theta*Omega;      %defined on page 61
+kappa=lambda*(siggma+(varphi+alppha)/(1-alppha));     %defined on page 63
+
+prior_mean=0.03;
+prior_std=0.02;
+log_prior_val=log(normpdf(kappa,prior_mean,prior_std));

examples/Occbin_example.mod

@@ -0,0 +1,139 @@
+/*
+ * This file shows how to solve an RBC model with two occasionally binding constraints:
+ *  1. The INEG constraint implements quadratic capital adjustment costs if investment
+ *      falls below its steady state. If investment is above steady state, there are no
+ *      adjustment costs
+ *  2. The IRR constraint implements irreversible investment. Investment cannot be lower
+ *      than a factor phi of its steady state.
+ *
+ * Notes:
+ *  - This mod-file is based on an example originally provided by Luca Guerrieri
+ *      and Matteo Iacoviello provided at https://www.matteoiacoviello.com/research_files/occbin_20140630.zip
+ *  - The INEG constraint should theoretically be log_Invest-log(steady_state(Invest))<0, but this will lead
+ *      to numerical issues. Instead we allow for a small negative value of <-0.000001
+ *
+ * Please note that the following copyright notice only applies to this Dynare
+ * implementation of the model.
+ */
+
+/*
+ * Copyright (C) 2021 Dynare Team
+ *
+ * This 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.
+ *
+ * It 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.
+ *
+ * For A copy of the GNU General Public License,
+ * see <https://www.gnu.org/licenses/>.
+ */
+
+
+var A           $A$         (long_name='TFP')
+    C           $C$         (long_name='consumption')
+    Invest      $I$         (long_name='investment')
+    K           $K$         (long_name='capital')
+    Lambda      $\lambda$   (long_name='Lagrange multiplier')
+    log_K       ${\hat K}$  (long_name='log capital')
+    log_Invest  ${\hat I}$  (long_name='log investment')
+    log_C       ${\hat C}$  (long_name='log consumption')
+    ;
+
+varexo epsilon $\varepsilon$        (long_name='TFP shock');
+
+parameters alpha    $\alpha$        (long_name='capital share')
+        delta       $\delta$        (long_name='depreciation')
+        beta        $\beta$         (long_name='discount factor')
+        sigma       $\sigma$        (long_name='risk aversion')
+        rho         $\rho$          (long_name='autocorrelation TFP')
+        phi         $\phi$          (long_name='irreversibility fraction of steady state investment')
+        psi         $\psi$          (long_name='capital adjustment cost')
+        ;
+
+beta=0.96;
+alpha=0.33;
+delta=0.10;
+sigma=2;
+rho = 0.9;
+phi = 0.975;
+psi = 5;
+
+
+model;
+// 1.
+[name='Euler', bind = 'INEG']
+-C^(-sigma)*(1+2*psi*(K/K(-1)-1)/K(-1))+ beta*C(+1)^(-sigma)*((1-delta)-2*psi*(K(+1)/K-1)*
+  (-K(+1)/K^2)+alpha*exp(A(+1))*K^(alpha-1))= -Lambda+beta*(1-delta)*Lambda(+1);
+
+[name='Euler', relax = 'INEG']
+-C^(-sigma) + beta*C(+1)^(-sigma)*(1-delta+alpha*exp(A(+1))*K^(alpha-1))= -Lambda+beta*(1-delta)*Lambda(+1);
+
+// 2.
+[name='Budget constraint',bind = 'INEG']
+C+K-(1-delta)*K(-1)+psi*(K/K(-1)-1)^2=exp(A)*K(-1)^(alpha);
+
+[name='Budget constraint',relax = 'INEG']
+C+K-(1-delta)*K(-1)=exp(A)*K(-1)^(alpha);
+
+// 3.
+[name='LOM capital']
+Invest = K-(1-delta)*K(-1);
+
+// 4.
+[name='investment',bind='IRR,INEG']
+(log_Invest - log(phi*steady_state(Invest))) = 0;
+[name='investment',relax='IRR']
+Lambda=0;
+[name='investment',bind='IRR',relax='INEG']
+(log_Invest - log(phi*steady_state(Invest))) = 0;
+
+// 5.
+[name='LOM TFP']
+A = rho*A(-1)+epsilon;
+
+// Definitions
+[name='Definition log capital']
+log_K=log(K);
+[name='Definition log consumption']
+log_C=log(C);
+[name='Definition log investment']
+log_Invest=log(Invest);
+end;
+
+occbin_constraints;
+name 'IRR'; bind log_Invest-log(steady_state(Invest))<log(phi); relax Lambda<0;
+name 'INEG'; bind log_Invest-log(steady_state(Invest))<-0.000001; %not exactly 0 for numerical reasons
+end;
+
+steady_state_model;
+K = ((1/beta-1+delta)/alpha)^(1/(alpha-1));
+C = -delta*K +K^alpha;
+Invest = delta*K;
+log_K = log(K);
+log_C = log(C);
+log_Invest = log(Invest);
+Lambda = 0;
+A=0;
+end;
+
+shocks;
+  var epsilon; stderr 0.015;
+end;
+
+steady;
+
+shocks(surprise);
+var epsilon;
+periods 1:9, 10, 50, 90, 130, 131:169;
+values -0.0001, -0.01,-0.02, 0.01, 0.02, 0;
+end;
+
+occbin_setup;
+occbin_solver(simul_periods=200,simul_check_ahead_periods=200);
+
+occbin_graph log_C epsilon Lambda log_K log_Invest A;

matlab/+occbin/DSGE_smoother.m

@@ -106,6 +106,7 @@ opts_simul.maxit = options_.occbin.smoother.maxit;
 opts_simul.waitbar = options_.occbin.smoother.waitbar;
 opts_simul.periods = options_.occbin.smoother.periods;
 opts_simul.check_ahead_periods = options_.occbin.smoother.check_ahead_periods;
+opts_simul.periodic_solution = options_.occbin.smoother.periodic_solution;
 opts_simul.full_output = options_.occbin.smoother.full_output;
 opts_simul.piecewise_only = options_.occbin.smoother.piecewise_only;
 % init_mode = options_.occbin.smoother.init_mode; % 0 = standard;  1 = unconditional frcsts zero shocks+smoothed states in each period

matlab/+occbin/set_default_options.m

@@ -73,14 +73,14 @@ if ismember(flag,{'likelihood','all'})
     options_occbin_.likelihood.init_binding_indicator = false(0);
     options_occbin_.likelihood.inversion_filter = false;
     options_occbin_.likelihood.IVF_shock_observable_mapping = [];
-    options_occbin_.likelihood.maxit = 50; % this is for occbin solver algo
+    options_occbin_.likelihood.maxit = 30; % this is for occbin solver algo
     options_occbin_.likelihood.max_number_of_iterations = 10; % this is for occbin_kalman_update loop
     options_occbin_.likelihood.periods = 100;
     options_occbin_.likelihood.check_ahead_periods=200;
-    options_occbin_.likelihood.periodic_solution=true;
+    options_occbin_.likelihood.periodic_solution=false;
     options_occbin_.likelihood.piecewise_only = true;
     options_occbin_.likelihood.restrict_state_space = true;
-    options_occbin_.likelihood.status=true;
+    options_occbin_.likelihood.status=true; %initialized to false in default_option_values
     options_occbin_.likelihood.use_updated_regime = true;
     options_occbin_.likelihood.waitbar=false;
 end
@@ -168,12 +168,11 @@ if ismember(flag,{'simul','all'})
     options_occbin_.simul.init_regime=[];
     options_occbin_.simul.init_binding_indicator=false(0);
     options_occbin_.simul.exo_pos=1:M_.exo_nbr;
-    options_occbin_.simul.local=true;
-    options_occbin_.simul.maxit=10;
+    options_occbin_.simul.maxit=30;
     options_occbin_.simul.max_periods=inf;
-    options_occbin_.simul.periods=30;
+    options_occbin_.simul.periods=100;
     options_occbin_.simul.check_ahead_periods=200;
-    options_occbin_.simul.periodic_solution=true;
+    options_occbin_.simul.periodic_solution=false;
     options_occbin_.simul.piecewise_only = false;
     options_occbin_.simul.reset_regime_in_new_period = false;
     options_occbin_.simul.restrict_state_space=false;
@@ -196,6 +195,7 @@ if ismember(flag,{'smoother','all'})
     options_occbin_.smoother.max_number_of_iterations = 10; % this is for smoother loop
     options_occbin_.smoother.periods = 100;
     options_occbin_.smoother.check_ahead_periods=200;
+    options_occbin_.smoother.periodic_solution=false;
     options_occbin_.smoother.piecewise_only = true;
     options_occbin_.smoother.plot = true;
     options_occbin_.smoother.status=true;

matlab/+occbin/setup.m

@@ -32,13 +32,23 @@ options_ = occbin.set_option(options_,options_occbin_,'simul.periods');
 options_ = occbin.set_option(options_,options_occbin_,'simul.curb_retrench');
 options_ = occbin.set_option(options_,options_occbin_,'simul.maxit');
 options_ = occbin.set_option(options_,options_occbin_,'simul.check_ahead_periods');
+options_ = occbin.set_option(options_,options_occbin_,'simul.debug');
+options_ = occbin.set_option(options_,options_occbin_,'simul.periodic_solution');
 options_ = occbin.set_option(options_,options_occbin_,'smoother.periods');
 options_ = occbin.set_option(options_,options_occbin_,'smoother.curb_retrench');
 options_ = occbin.set_option(options_,options_occbin_,'smoother.maxit');
 options_ = occbin.set_option(options_,options_occbin_,'smoother.check_ahead_periods');
+options_ = occbin.set_option(options_,options_occbin_,'smoother.debug');
+options_ = occbin.set_option(options_,options_occbin_,'smoother.periodic_solution');
+options_ = occbin.set_option(options_,options_occbin_,'smoother.inversion_filter');
+options_ = occbin.set_option(options_,options_occbin_,'likelihood.periods');
+options_ = occbin.set_option(options_,options_occbin_,'likelihood.curb_retrench');
+options_ = occbin.set_option(options_,options_occbin_,'likelihood.maxit');
+options_ = occbin.set_option(options_,options_occbin_,'likelihood.check_ahead_periods');
+options_ = occbin.set_option(options_,options_occbin_,'likelihood.periodic_solution');
+options_ = occbin.set_option(options_,options_occbin_,'likelihood.max_number_of_iterations');
 options_ = occbin.set_option(options_,options_occbin_,'filter.use_relaxation');
 options_ = occbin.set_option(options_,options_occbin_,'likelihood.inversion_filter');
-options_ = occbin.set_option(options_,options_occbin_,'smoother.inversion_filter');
 
 if isfield(M_,'surprise_shocks') && ~isempty(M_.surprise_shocks)
     temp=zeros(max(cat(2,M_.surprise_shocks.periods)),M_.exo_nbr);

matlab/+occbin/solver.m

@@ -77,7 +77,6 @@ if ~options_.occbin.simul.piecewise_only
     out.linear    = out.linear + out.ys';
 end
 out.piecewise = out.piecewise+ out.ys';
-out.exo_simul = options_.occbin.simul.SHOCKS;
 out.exo_pos = options_.occbin.simul.exo_pos;
 
 oo_.occbin=out;

matlab/GetPosteriorParametersStatistics.m

@@ -63,7 +63,7 @@ clear record;
 header_width = row_header_width(M_, estim_params_, bayestopt_);
 hpd_interval=[num2str(options_.mh_conf_sig*100), '% HPD interval'];
 tit2 = sprintf('%-*s %12s %12s %23s %8s %12s\n',header_width,' ','prior mean','post. mean',hpd_interval,'prior','pstdev');
-pformat =      '%-*s %12.3f % 12.4f %11.4f %11.4f %7s %12.4f';
+pformat =      '%-*s %12.3f % 12.4f %11.4f %11.4f %8s %12.4f';
 
 skipline(2)
 disp('ESTIMATION RESULTS')

matlab/check_particle_filter_options.m

@@ -111,12 +111,6 @@ if ~isempty(particle_options.particle_filter_options)
                 else
                     particle_options.liu_west_delta=options_list{i,2};
                 end
-            case 'liu_west_chol_sigma_bar'
-                if options_list{i,2} <= 0
-                    error('check_particle_filter_options:: the liu_west_chol_sigma_bar option takes a positive argument');
-                else
-                    particle_options.liu_west_chol_sigma_bar=options_list{i,2};
-                end
                 
             otherwise
                 warning(['check_particle_filter_options: Unknown option (' options_list{i,1} ')!'])

matlab/default_option_values.m

@@ -286,7 +286,6 @@ particle.mixture_structural_shocks = 1 ;
 particle.mixture_measurement_shocks = 1 ;
 % Online approach
 particle.liu_west_delta = 0.99 ;
-particle.liu_west_chol_sigma_bar = .01 ;
 % Options for setting the weights in conditional particle filters.
 particle.cpf_weights_method.amisanotristani = true;
 particle.cpf_weights_method.murrayjonesparslow = false;

matlab/display_estimation_results_table.m

@@ -51,10 +51,10 @@ tstath = abs(xparam1)./stdh;
 
 header_width = row_header_width(M_, estim_params_, bayestopt_);
 if strcmp(field_name,'posterior')
-    tit1 = sprintf('%-*s %7s %8s %7s %4s %6s\n', header_width-2, ' ', 'prior mean', ...
+    tit1 = sprintf('%-*s %10s %8s %7s %6s %6s\n', header_width, ' ', 'prior mean', ...
                    'mode', 's.d.', 'prior', 'pstdev');
 else
-    tit1 = sprintf('%-*s %10s %7s %6s\n', header_width-2, ' ', 'Estimate', 's.d.', 't-stat');
+    tit1 = sprintf('%-*s %10s %7s %6s\n', header_width, ' ', 'Estimate', 's.d.', 't-stat');
 end
 if np
     ip = nvx+nvn+ncx+ncn+1;
@@ -63,13 +63,13 @@ if np
     for i=1:np
         name = bayestopt_.name{ip};
         if strcmp(field_name,'posterior')
-            fprintf('%-*s %7.3f %8.4f %7.4f %4s %6.4f \n', ...
+            fprintf('%-*s %10.4f %8.4f %7.4f %6s %6.4f \n', ...
                     header_width,name, ...
                     bayestopt_.p1(ip),xparam1(ip),stdh(ip), ...
                     pnames{bayestopt_.pshape(ip)+1}, ...
                     bayestopt_.p2(ip));
         else
-            fprintf('%-*s %8.4f %7.4f %7.4f \n', ...
+            fprintf('%-*s %10.4f %7.4f %7.4f \n', ...
                     header_width, name, xparam1(ip), stdh(ip), tstath(ip));
         end
         eval(['oo_.' field_name '_mode.parameters.' name ' = xparam1(ip);']);
@@ -86,12 +86,12 @@ if nvx
         k = estim_params_.var_exo(i,1);
         name = M_.exo_names{k};
         if strcmp(field_name,'posterior')
-            fprintf('%-*s %7.3f %8.4f %7.4f %4s %6.4f \n', ...
+            fprintf('%-*s %10.4f %8.4f %7.4f %6s %6.4f \n', ...
                     header_width, name, bayestopt_.p1(ip), xparam1(ip), ...
                     stdh(ip), pnames{bayestopt_.pshape(ip)+1}, ...
                     bayestopt_.p2(ip));
         else
-            fprintf('%-*s %8.4f %7.4f %7.4f \n', header_width, name, xparam1(ip), stdh(ip), tstath(ip));
+            fprintf('%-*s %10.4f %7.4f %7.4f \n', header_width, name, xparam1(ip), stdh(ip), tstath(ip));
         end
         M_.Sigma_e(k,k) = xparam1(ip)*xparam1(ip);
         eval(['oo_.' field_name '_mode.shocks_std.' name ' = xparam1(ip);']);
@@ -107,13 +107,13 @@ if nvn
     for i=1:nvn
         name = options_.varobs{estim_params_.nvn_observable_correspondence(i,1)};
         if strcmp(field_name,'posterior')
-            fprintf('%-*s %7.3f %8.4f %7.4f %4s %6.4f \n', ...
+            fprintf('%-*s %10.4f %8.4f %7.4f %6s %6.4f \n', ...
                     header_width, name, bayestopt_.p1(ip), ...
                     xparam1(ip), stdh(ip), ...
                     pnames{bayestopt_.pshape(ip)+1}, ...
                     bayestopt_.p2(ip));
         else
-            fprintf('%-*s %8.4f %7.4f %7.4f \n', header_width, name, xparam1(ip), ...
+            fprintf('%-*s %10.4f %7.4f %7.4f \n', header_width, name, xparam1(ip), ...
                     stdh(ip), tstath(ip))
         end
         eval(['oo_.' field_name '_mode.measurement_errors_std.' name ' = xparam1(ip);']);
@@ -133,11 +133,11 @@ if ncx
         name = sprintf('%s,%s', M_.exo_names{k1}, M_.exo_names{k2});
         NAME = sprintf('%s_%s', M_.exo_names{k1}, M_.exo_names{k2});
         if strcmp(field_name, 'posterior')
-            fprintf('%-*s %7.3f %8.4f %7.4f %4s %6.4f \n', ...
+            fprintf('%-*s %10.4f %8.4f %7.4f %6s %6.4f \n', ...
                     header_width, name, bayestopt_.p1(ip), xparam1(ip), stdh(ip),  ...
                     pnames{bayestopt_.pshape(ip)+1}, bayestopt_.p2(ip));
         else
-            fprintf('%-*s %8.4f %7.4f %7.4f \n', header_width,name, xparam1(ip), ...
+            fprintf('%-*s %10.4f %7.4f %7.4f \n', header_width,name, xparam1(ip), ...
                     stdh(ip), tstath(ip));
         end
         M_.Sigma_e(k1,k2) = xparam1(ip)*sqrt(M_.Sigma_e(k1,k1)*M_.Sigma_e(k2,k2));
@@ -159,11 +159,11 @@ if ncn
         name = sprintf('%s,%s', M_.endo_names{k1}, M_.endo_names{k2});
         NAME = sprintf('%s_%s', M_.endo_names{k1}, M_.endo_names{k2});
         if strcmp(field_name,'posterior')
-            fprintf('%-*s %7.3f %8.4f %7.4f %4s %6.4f \n', ...
+            fprintf('%-*s %10.4f %8.4f %7.4f %6s %6.4f \n', ...
                     header_width, name, bayestopt_.p1(ip), xparam1(ip), stdh(ip), ...
                     pnames{bayestopt_.pshape(ip)+1}, bayestopt_.p2(ip));
         else
-            fprintf('%-*s %8.4f %7.4f %7.4f \n',header_width, name, xparam1(ip), ...
+            fprintf('%-*s %10.4f %7.4f %7.4f \n',header_width, name, xparam1(ip), ...
                     stdh(ip), tstath(ip));
         end
         eval(['oo_.' field_name '_mode.measurement_errors_corr.' NAME ' = xparam1(ip);']);

matlab/dr_block.m

@@ -61,11 +61,18 @@ end
 %verbose = options_.verbosity;
 if options_.order > 1
     error('2nd and 3rd order approximation not implemented with block option')
+else
+    if options_.loglinear
+        error('The loglinear option is not yet supported in first order approximation for a block decomposed model');
+    end
 end
 
 z = repmat(dr.ys,1,M_.maximum_lead + M_.maximum_lag + 1);
 zx = repmat([oo_.exo_simul oo_.exo_det_simul],M_.maximum_lead + M_.maximum_lag + 1, 1);
-
+if isempty(zx)
+    zx = [repmat(oo_.exo_steady_state',M_.maximum_lead + M_.maximum_lag + 1,1) repmat(oo_.exo_det_steady_state',M_.maximum_lead + M_.maximum_lag + 1,1)];
+end
+    
 if ~isfield(M_,'block_structure')
     error('Option ''block'' has not been specified')
 end

matlab/dsge_likelihood.m

@@ -670,7 +670,7 @@ singularity_has_been_detected = false;
 if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter
     if no_missing_data_flag && ~DynareOptions.occbin.likelihood.status
         if DynareOptions.block
-            LIK = block_kalman_filter(T,R,Q,H,Pstar,Y,start,Z,kalman_tol,riccati_tol, Model.nz_state_var, Model.n_diag, Model.nobs_non_statevar);
+            [LIK,lik] = block_kalman_filter(T,R,Q,H,Pstar,Y,start,Z,kalman_tol,riccati_tol, Model.nz_state_var, Model.n_diag, Model.nobs_non_statevar);
         elseif DynareOptions.fast_kalman_filter
             if diffuse_periods
                 %kalman_algo==3 requires no diffuse periods (stationary

matlab/evaluate_planner_objective.m

@@ -6,8 +6,8 @@ function planner_objective_value = evaluate_planner_objective(M_,options_,oo_)
 %   oo_:       (structure) output results
 % OUTPUT
 %  planner_objective_value (double)
-% 
-%Returns a vector containing first order or second-order approximations of 
+%
+%Returns a vector containing first order or second-order approximations of
 % - the unconditional expectation of the planner's objective function
 % - the conditional expectation of the planner's objective function starting from the non-stochastic steady state and allowing for future shocks
 % depending on the value of options_.order.
@@ -52,7 +52,7 @@ function planner_objective_value = evaluate_planner_objective(M_,options_,oo_)
 % W(y,0,1) = Wbar + 0.5*Wss
 
 % In the discretionary case, the model is assumed to be linear and the utility is assumed to be linear-quadratic. This changes 2 aspects of the results delinated above:
-% 1) the second-order derivatives of the policy and transition functions h and g are zero. 
+% 1) the second-order derivatives of the policy and transition functions h and g are zero.
 % 2) the unconditional expectation of states coincides with its steady-state, which entails E(yhat) = 0
 % Therefore, the unconditional welfare can now be approximated as
 % E(W) = (1 - beta)^{-1} ( Ubar + 0.5 ( U_xx h_y^2 E(yhat^2) + U_xx h_u^2 E(u^2) )
@@ -91,47 +91,47 @@ planner_objective_value = zeros(2,1);
 if options_.ramsey_policy
     if oo_.gui.ran_perfect_foresight
         T = size(oo_.endo_simul,2);
-        [U_term] = feval([M_.fname '.objective.static'],oo_.endo_simul(:,T),oo_.exo_simul(T,:), M_.params);
+        [U_term] = feval([M_.fname '.objective.static'],oo_.endo_simul(:,T-M_.maximum_lead),oo_.exo_simul(T-M_.maximum_lead,:), M_.params);
         EW = U_term/(1-beta);
         W = EW;
-        for t=T:-1:2
+        for t=T-M_.maximum_lead:-1:1+M_.maximum_lag
             [U] = feval([M_.fname '.objective.static'],oo_.endo_simul(:,t),oo_.exo_simul(t,:), M_.params);
             W = U + beta*W;
         end
         planner_objective_value(1) = EW;
         planner_objective_value(2) = W;
     else
-    ys = oo_.dr.ys;
+        ys = oo_.dr.ys;
         if options_.order == 1 || M_.hessian_eq_zero
             [U] = feval([M_.fname '.objective.static'],ys,zeros(1,exo_nbr), M_.params);
             planner_objective_value(1) = U/(1-beta);
-            planner_objective_value(2) = U/(1-beta);  
+            planner_objective_value(2) = U/(1-beta);
         elseif options_.order == 2 && ~M_.hessian_eq_zero
             [U,Uy,Uyy] = feval([M_.fname '.objective.static'],ys,zeros(1,exo_nbr), M_.params);
-
+            
             Gy = dr.ghx(nstatic+(1:nspred),:);
             Gu = dr.ghu(nstatic+(1:nspred),:);
             Gyy = dr.ghxx(nstatic+(1:nspred),:);
             Gyu = dr.ghxu(nstatic+(1:nspred),:);
             Guu = dr.ghuu(nstatic+(1:nspred),:);
             Gss = dr.ghs2(nstatic+(1:nspred),:);
-
+            
             gy(dr.order_var,:) = dr.ghx;
             gu(dr.order_var,:) = dr.ghu;
             gyy(dr.order_var,:) = dr.ghxx;
             gyu(dr.order_var,:) = dr.ghxu;
             guu(dr.order_var,:) = dr.ghuu;
             gss(dr.order_var,:) = dr.ghs2;
-
+            
             Uyy = full(Uyy);
-
+            
             Uyygygy = A_times_B_kronecker_C(Uyy,gy,gy);
             Uyygugu = A_times_B_kronecker_C(Uyy,gu,gu);
-
+            
             %% Unconditional welfare
-
+            
             old_noprint = options_.noprint;
-
+            
             if ~old_noprint
                 options_.noprint = 1;
             end
@@ -143,25 +143,25 @@ if options_.ramsey_policy
             if ~old_noprint
                 options_.noprint = 0;
             end
-
+            
             oo_.mean(isnan(oo_.mean)) = options_.huge_number;
             oo_.var(isnan(oo_.var)) = options_.huge_number;
-
+            
             Ey = oo_.mean;
             Eyhat = Ey - ys(dr.order_var(nstatic+(1:nspred)));
-
+            
             var_corr = Eyhat*Eyhat';
             Eyhatyhat = oo_.var(:) + var_corr(:);
             Euu = M_.Sigma_e(:);
-
+            
             EU = U + Uy*gy*Eyhat + 0.5*((Uyygygy + Uy*gyy)*Eyhatyhat + (Uyygugu + Uy*guu)*Euu + Uy*gss);
             EW = EU/(1-beta);
-
+            
             %% Conditional welfare starting from the non-stochastic steady-state
-
+            
             Wbar = U/(1-beta);
             Wy = Uy*gy/(eye(nspred)-beta*Gy);
-
+            
             if isempty(options_.qz_criterium)
                 options_.qz_criterium = 1+1e-6;
             end
@@ -171,7 +171,7 @@ if options_.ramsey_policy
             Wuu = Uyygugu + Uy*guu + beta*(Wyygugu + Wy*Guu);
             Wss = (Uy*gss + beta*(Wy*Gss + Wuu*M_.Sigma_e(:)))/(1-beta);
             W = Wbar + 0.5*Wss;
-
+            
             planner_objective_value(1) = EW;
             planner_objective_value(2) = W;
         else
@@ -184,23 +184,23 @@ if options_.ramsey_policy
     end
 elseif options_.discretionary_policy
     ys = oo_.dr.ys;
-        
+    
     [U,Uy,Uyy] = feval([M_.fname '.objective.static'],ys,zeros(1,exo_nbr), M_.params);
-
+    
     Gy = dr.ghx(nstatic+(1:nspred),:);
     Gu = dr.ghu(nstatic+(1:nspred),:);
     gy(dr.order_var,:) = dr.ghx;
     gu(dr.order_var,:) = dr.ghu;
-
+    
     Uyy = full(Uyy);
-
+    
     Uyygygy = A_times_B_kronecker_C(Uyy,gy,gy);
     Uyygugu = A_times_B_kronecker_C(Uyy,gu,gu);
-
+    
     %% Unconditional welfare
-
+    
     old_noprint = options_.noprint;
-
+    
     if ~old_noprint
         options_.noprint = 1;
     end
@@ -209,36 +209,33 @@ elseif options_.discretionary_policy
     if ~old_noprint
         options_.noprint = 0;
     end
-
+    
     oo_.mean(isnan(oo_.mean)) = options_.huge_number;
     oo_.var(isnan(oo_.var)) = options_.huge_number;
-
+    
     Ey = oo_.mean;
     Eyhat = Ey - ys(dr.order_var(nstatic+(1:nspred)));
-
+    
     var_corr = Eyhat*Eyhat';
     Eyhatyhat = oo_.var(:) + var_corr(:);
     Euu = M_.Sigma_e(:);
-
+    
     EU = U + Uy*gy*Eyhat + 0.5*(Uyygygy*Eyhatyhat + Uyygugu*Euu);
     EW = EU/(1-beta);
-
-
+    
+    
     %% Conditional welfare starting from the non-stochastic steady-state
-
+    
     Wbar = U/(1-beta);
     Wy = Uy*gy/(eye(nspred)-beta*Gy);
-
-    if isempty(options_.qz_criterium)
-        options_.qz_criterium = 1+1e-6;
-    end
+    
     %solve Lyapunuv equation Wyy=gy'*Uyy*gy+beta*Gy'Wyy*Gy
     Wyy = reshape(lyapunov_symm(sqrt(beta)*Gy',reshape(Uyygygy,nspred,nspred),options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold, 3, options_.debug),1,nspred*nspred);
     Wyygugu = A_times_B_kronecker_C(Wyy,Gu,Gu);
     Wuu = Uyygugu + beta*Wyygugu;
     Wss = beta*Wuu*M_.Sigma_e(:)/(1-beta);
     W = Wbar + 0.5*Wss;
-
+    
     planner_objective_value(1) = EW;
     planner_objective_value(2) = W;
 end
@@ -250,7 +247,7 @@ if ~options_.noprint
         else
             fprintf('\nApproximated value of unconditional welfare:  %10.8f\n', planner_objective_value(1))
             fprintf('\nApproximated value of conditional welfare:  %10.8f\n', planner_objective_value(2))
-        end 
+        end
     elseif options_.discretionary_policy
         fprintf('\nApproximated value of unconditional welfare with discretionary policy:  %10.8f\n', planner_objective_value(1))
         fprintf('\nApproximated value of conditional welfare with discretionary policy:  %10.8f\n', planner_objective_value(2))

matlab/execute_prior_posterior_function.m

@@ -68,12 +68,9 @@ elseif strcmpi(type,'prior')
         else
             error('The prior distributions are not properly set up.')
         end
-        if exist([M_.fname '_prior_restrictions.m'])
-            error('prior_function currently does not support endogenous prior restrictions.')
-            %options_.plot_priors=0;
-            %[~,~,~,~, M_, options_, oo_, ~, ~] = ...
-            %    dynare_estimation_init(M_.endo_names, M_.fname, 1, M_, options_, oo_, estim_params_, bayestopt_);
-        end
+    end
+    if exist([M_.fname '_prior_restrictions.m'])
+        warning('prior_function currently does not support endogenous prior restrictions. They will be ignored. Consider using a prior_function with nobs=1.')
     end
     prior_draw(bayestopt_, options_.prior_trunc);
 else

matlab/method_of_moments/method_of_moments.m

@@ -91,8 +91,8 @@ function [oo_, options_mom_, M_] = method_of_moments(bayestopt_, options_, oo_,
 % - [ ] add option to use autocorrelations (we have useautocorr in identification toolbox already)
 % - [ ] SMM with extended path
 % - [ ] deal with measurement errors (once @wmutschl has implemented this in identification toolbox)
-% - [ ] improve check for duplicate moments by using the cellfun and unique functions
 % - [ ] dirname option to save output to different directory not yet implemented
+% - [ ] display scaled moments
 
 % The TeX option crashes MATLAB R2014a run with "-nodisplay" option
 % (as is done from the testsuite).
@@ -152,7 +152,7 @@ if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_meth
     options_mom_.mom = set_default_option(options_mom_.mom,'bartlett_kernel_lag',20);               % bandwith in optimal weighting matrix
     options_mom_.mom = set_default_option(options_mom_.mom,'penalized_estimator',false);            % include deviation from prior mean as additional moment restriction and use prior precision as weight
     options_mom_.mom = set_default_option(options_mom_.mom,'verbose',false);                        % display and store intermediate estimation results
-    options_mom_.mom = set_default_option(options_mom_.mom,'weighting_matrix',{'DIAGONAL'; 'DIAGONAL'});   % weighting matrix in moments distance objective function at each iteration of estimation;
+    options_mom_.mom = set_default_option(options_mom_.mom,'weighting_matrix',{'DIAGONAL'; 'OPTIMAL'});   % weighting matrix in moments distance objective function at each iteration of estimation;
                                                                                                            % possible values are 'OPTIMAL', 'IDENTITY_MATRIX' ,'DIAGONAL' or a filename. Size of cell determines stages in iterated estimation.
     options_mom_.mom = set_default_option(options_mom_.mom,'weighting_matrix_scaling_factor',1);    % scaling of weighting matrix in objective function
     options_mom_.mom = set_default_option(options_mom_.mom,'se_tolx',1e-5);                         % step size for numerical computation of standard errors
@@ -167,10 +167,10 @@ if strcmp(options_mom_.mom.mom_method,'SMM')
     options_mom_.mom = set_default_option(options_mom_.mom,'burnin',500);                           % number of periods dropped at beginning of simulation
     options_mom_.mom = set_default_option(options_mom_.mom,'bounded_shock_support',false);          % trim shocks in simulation to +- 2 stdev
     options_mom_.mom = set_default_option(options_mom_.mom,'seed',24051986);                        % seed used in simulations
-    options_mom_.mom = set_default_option(options_mom_.mom,'simulation_multiple',5);                % multiple of the data length used for simulation
+    options_mom_.mom = set_default_option(options_mom_.mom,'simulation_multiple',7);                % multiple of the data length used for simulation
     if options_mom_.mom.simulation_multiple < 1
         fprintf('The simulation horizon is shorter than the data. Dynare resets the simulation_multiple to 5.\n')
-        options_mom_.mom.simulation_multiple = 5;
+        options_mom_.mom.simulation_multiple = 7;
     end
 end
 if strcmp(options_mom_.mom.mom_method,'GMM')
@@ -215,7 +215,11 @@ end
 
 % Optimization options that can be set by the user in the mod file, otherwise default values are provided
 options_mom_ = set_default_option(options_mom_,'huge_number',1e7);               % value for replacing the infinite bounds on parameters by finite numbers. Used by some optimizers for numerical reasons
-options_mom_ = set_default_option(options_mom_,'mode_compute',13);               % specifies the optimizer for minimization of moments distance
+if (isoctave && user_has_octave_forge_package('optim')) || (~isoctave && user_has_matlab_license('optimization_toolbox'))
+    options_mom_ = set_default_option(options_mom_,'mode_compute',13);               % specifies lsqnonlin as default optimizer for minimization of moments distance
+else
+    options_mom_ = set_default_option(options_mom_,'mode_compute',4);               % specifies csminwel as fallback default option for minimization of moments distance
+end
 options_mom_ = set_default_option(options_mom_,'additional_optimizer_steps',[]); % vector of additional mode-finders run after mode_compute
 options_mom_ = set_default_option(options_mom_,'optim_opt',[]);                  % a list of NAME and VALUE pairs to set options for the optimization routines. Available options depend on mode_compute
 options_mom_ = set_default_option(options_mom_,'silent_optimizer',false);        % run minimization of moments distance silently without displaying results or saving files in between
@@ -384,85 +388,66 @@ for i=1:options_mom_.obs_nbr
 end
 
 % -------------------------------------------------------------------------
-% Step 2: Checks and transformations for matched moments structure (preliminary)
+% Step 2: Checks and transformations for matched_moments structure
 % -------------------------------------------------------------------------
+M_.matched_moments_orig = M_.matched_moments; %save original structure
 
-% Initialize indices
-options_mom_.mom.index.E_y       = false(options_mom_.obs_nbr,1);                      %unconditional first order product moments
-options_mom_.mom.index.E_yy      = false(options_mom_.obs_nbr,options_mom_.obs_nbr);   %unconditional second order product moments
-options_mom_.mom.index.E_yyt     = false(options_mom_.obs_nbr,options_mom_.obs_nbr,0); %unconditional temporal second order product moments
-options_mom_.mom.index.E_y_pos   = zeros(options_mom_.obs_nbr,1);                      %position in matched moments block
-options_mom_.mom.index.E_yy_pos  = zeros(options_mom_.obs_nbr,options_mom_.obs_nbr);   %position in matched moments block
-options_mom_.mom.index.E_yyt_pos = zeros(options_mom_.obs_nbr,options_mom_.obs_nbr,0); %position in matched moments block
+% higher-order product moments not supported yet for GMM
+if strcmp(options_mom_.mom.mom_method, 'GMM') && any(cellfun(@sum,M_.matched_moments(:,3))> 2)
+    error('method_of_moments: GMM does not yet support product moments higher than 2. Change row %d in ''matched_moments'' block.',jm);
+end
 
+% check for duplicate moment conditions
 for jm=1:size(M_.matched_moments,1)
-    % higher-order product moments not supported yet for GMM
-    if strcmp(options_mom_.mom.mom_method, 'GMM') && sum(M_.matched_moments{jm,3}) > 2
-        error('method_of_moments: GMM does not yet support product moments higher than 2. Change row %d in ''matched_moments'' block.',jm);
-    end
-    % Check if declared variables are also observed (needed as otherwise the dataset variables won't coincide)
-    if any(~ismember(oo_.dr.inv_order_var(M_.matched_moments{jm,1})', oo_.dr.obs_var))
-        error('method_of_moments: Variables in row %d in ''matched_moments'' block need to be declared as VAROBS.', jm)
-    end
-    
-    if strcmp(options_mom_.mom.mom_method, 'GMM')
-    % Check (for now) that only lags are declared
-        if any(M_.matched_moments{jm,2}>0)
-            error('method_of_moments: Leads in row %d in the ''matched_moments'' block are not supported for GMM, shift the moments and declare only lags.', jm)
-        end
-        % Check (for now) that first declared variable has zero lag
-        if M_.matched_moments{jm,2}(1)~=0
-            error('method_of_moments: The first variable declared in row %d in the ''matched_moments'' block is not allowed to have a lead or lag for GMM;\n                   reorder the variables in the row such that the first variable has zero lag!',jm)
-        end
-    end
-    vars = oo_.dr.inv_order_var(M_.matched_moments{jm,1})';
-    if sum(M_.matched_moments{jm,3}) == 1
-        % First-order product moment
-        vpos = (oo_.dr.obs_var == vars);
-        options_mom_.mom.index.E_y(vpos,1) = true;
-        options_mom_.mom.index.E_y_pos(vpos,1) = jm;
-        M_.matched_moments{jm,4}=['E(',M_.endo_names{M_.matched_moments{jm,1}},')'];
-        M_.matched_moments{jm,5}=['$E(',M_.endo_names_tex{M_.matched_moments{jm,1}},')$'];
-    elseif sum(M_.matched_moments{jm,3}) == 2
-        % Second-order product moment
-        idx1 = (oo_.dr.obs_var == vars(1));
-        idx2 = (oo_.dr.obs_var == vars(2));
-        lag1 = M_.matched_moments{jm,2}(1);
-        lag2 = M_.matched_moments{jm,2}(2);
-        if lag1==0 && lag2==0 % contemporaneous covariance matrix
-            options_mom_.mom.index.E_yy(idx1,idx2) = true;
-            options_mom_.mom.index.E_yy(idx2,idx1) = true;
-            options_mom_.mom.index.E_yy_pos(idx1,idx2) = jm;
-            options_mom_.mom.index.E_yy_pos(idx2,idx1) = jm;
-            M_.matched_moments{jm,4}=['E(',M_.endo_names{M_.matched_moments{jm,1}(1)},',',M_.endo_names{M_.matched_moments{jm,1}(2)},')'];
-            M_.matched_moments{jm,5}=['$E({',M_.endo_names_tex{M_.matched_moments{jm,1}(1)},'}_t,{',M_.endo_names_tex{M_.matched_moments{jm,1}(1)},'}_t)$'];
-        elseif lag1==0 && lag2 < 0
-            options_mom_.mom.index.E_yyt(idx1,idx2,-lag2) = true;
-            options_mom_.mom.index.E_yyt_pos(idx1,idx2,-lag2) = jm;
-            M_.matched_moments{jm,4}=['E(',M_.endo_names{M_.matched_moments{jm,1}(1)},',',M_.endo_names{M_.matched_moments{jm,1}(2)},'(',num2str(lag2),'))'];
-            M_.matched_moments{jm,5}=['$E({',M_.endo_names_tex{M_.matched_moments{jm,1}(1)},'}_t\times{',M_.endo_names_tex{M_.matched_moments{jm,1}(1)},'_{t',num2str(lag2) ,'})$'];
+    % expand powers to vector of ones
+    if any(M_.matched_moments{jm,3}>1)
+        tmp1=[]; tmp2=[]; tmp3=[];
+        for jjm=1:length(M_.matched_moments{jm,3})
+            tmp1 = [tmp1 repmat(M_.matched_moments{jm,1}(jjm),[1 M_.matched_moments{jm,3}(jjm)]) ];
+            tmp2 = [tmp2 repmat(M_.matched_moments{jm,2}(jjm),[1 M_.matched_moments{jm,3}(jjm)]) ];
+            tmp3 = [tmp3 repmat(1,[1 M_.matched_moments{jm,3}(jjm)]) ];
         end
+        M_.matched_moments{jm,1} = tmp1;
+        M_.matched_moments{jm,2} = tmp2;
+        M_.matched_moments{jm,3} = tmp3;
     end
+    % shift time structure to focus only on lags
+    M_.matched_moments{jm,2} = M_.matched_moments{jm,2} - max(M_.matched_moments{jm,2});
+    % sort such that t=0 variable comes first
+    [M_.matched_moments{jm,2},idx_sort] = sort(M_.matched_moments{jm,2},'descend');
+    M_.matched_moments{jm,1} = M_.matched_moments{jm,1}(idx_sort);
+    M_.matched_moments{jm,3} = M_.matched_moments{jm,3}(idx_sort);
 end
 
-%Remove duplicate elements
-UniqueMomIdx = [nonzeros(options_mom_.mom.index.E_y_pos); nonzeros(tril(options_mom_.mom.index.E_yy_pos)); nonzeros(options_mom_.mom.index.E_yyt_pos)];
-DuplicateMoms = setdiff(1:size(M_.matched_moments,1),UniqueMomIdx);
-if ~isempty(DuplicateMoms)
-    fprintf('Found and removed duplicate declared moments in ''matched_moments'' block in rows: %s.\n',num2str(DuplicateMoms))
+% find duplicate rows in cell array by making groups according to powers as we can then use cell2mat for the unique function
+powers = cellfun(@sum,M_.matched_moments(:,3))';
+UniqueMomIdx = [];
+for jpow = unique(powers)
+    idx1 = find(powers==jpow);
+    [~,idx2] = unique(cell2mat(M_.matched_moments(idx1,:)),'rows');
+    UniqueMomIdx = [UniqueMomIdx idx1(idx2)];
 end
-%reorder M_.matched_moments to be compatible with options_mom_.mom.index
-M_.matched_moments = M_.matched_moments(UniqueMomIdx,:);
-if strcmp(options_mom_.mom.mom_method,'SMM')
-    options_mom_.mom=rmfield(options_mom_.mom,'index');
+
+% remove duplicate elements
+DuplicateMoms = setdiff(1:size(M_.matched_moments_orig,1),UniqueMomIdx);
+if ~isempty(DuplicateMoms)
+    fprintf('Found and removed duplicate declared moments in ''matched_moments'' block in rows:\n %s.\n',num2str(DuplicateMoms))
+    fprintf('Dynare will continue with remaining moment conditions\n');
+end
+
+if strcmp(options_mom_.mom.mom_method, 'SMM')
+    % for SMM we can keep the original structure but get rid of duplicate moments
+    M_.matched_moments = M_.matched_moments_orig(sort(UniqueMomIdx),:);
+elseif strcmp(options_mom_.mom.mom_method, 'GMM')
+    % for GMM we use the transformed matched_moments structure
+    M_.matched_moments = M_.matched_moments(sort(UniqueMomIdx),:);
 end
 
 % Check if both prefilter and first moments were specified
-options_mom_.mom.first_moment_indicator = find(cellfun(@(x) sum(abs(x))==1,M_.matched_moments(:,3)))';
-if options_mom_.prefilter && ~isempty(options_mom_.mom.first_moment_indicator)
-    fprintf('Centered moments requested (prefilter option is set); therefore, ignore declared first moments in ''matched_moments'' block in rows: %u.\n',options_mom_.mom.first_moment_indicator');
-    M_.matched_moments(options_mom_.mom.first_moment_indicator,:)=[]; %remove first moments entries
-    options_mom_.mom.first_moment_indicator = [];
+first_moment_indicator = find(cellfun(@(x) sum(abs(x))==1,M_.matched_moments(:,3)))';
+if options_mom_.prefilter && ~isempty(first_moment_indicator)
+    fprintf('Centered moments requested (prefilter option is set); therefore, ignore declared first moments in ''matched_moments'' block.\n');
+    M_.matched_moments(first_moment_indicator,:)=[]; %remove first moments entries
 end
 options_mom_.mom.mom_nbr = size(M_.matched_moments,1);
 
@@ -497,7 +482,7 @@ end
 bayestopt_laplace=bayestopt_;
 
 % Check on specified priors and penalized estimation
-if any(bayestopt_laplace.pshape > 0) % prior specified, not ML
+if any(bayestopt_laplace.pshape > 0) % prior specified
     if ~options_mom_.mom.penalized_estimator
         fprintf('\nPriors were specified, but the penalized_estimator-option was not set.\n')
         fprintf('Dynare sets penalized_estimator to 1. Conducting %s with penalty.\n',options_mom_.mom.mom_method)
@@ -855,7 +840,6 @@ if size(options_mom_.mom.weighting_matrix,1)>1 && ~(any(strcmpi('diagonal',optio
     fprintf('\nYou did not specify the use of an optimal or diagonal weighting matrix. There is no point in running an iterated method of moments.\n')
 end
 
-optim_opt0 = options_mom_.optim_opt; % store original options set by user
 for stage_iter=1:size(options_mom_.mom.weighting_matrix,1)
     fprintf('Estimation stage %u\n',stage_iter);
     Woptflag = false;
@@ -972,15 +956,39 @@ end
 % -------------------------------------------------------------------------
 % Step 9: Display estimation results
 % -------------------------------------------------------------------------
-title = ['Data moments and model moments (',options_mom_.mom.mom_method,')'];
-headers = {'Moment','Data','Model','% dev. target'};
-labels= M_.matched_moments(:,4);
-data_mat=[oo_.mom.data_moments oo_.mom.model_moments 100*abs((oo_.mom.model_moments-oo_.mom.data_moments)./oo_.mom.data_moments)];
+title = ['Comparison of data moments and model moments (',options_mom_.mom.mom_method,')'];
+headers = {'Moment','Data','Model'};
+for jm = 1:size(M_.matched_moments,1)
+    lables_tmp = 'E[';
+    lables_tmp_tex = 'E \left[ ';
+    for jvar = 1:length(M_.matched_moments{jm,1})
+        lables_tmp = [lables_tmp M_.endo_names{M_.matched_moments{jm,1}(jvar)}];
+        lables_tmp_tex = [lables_tmp_tex, '{', M_.endo_names_tex{M_.matched_moments{jm,1}(jvar)}, '}'];
+        if M_.matched_moments{jm,2}(jvar) ~= 0
+            lables_tmp = [lables_tmp, '(', num2str(M_.matched_moments{jm,2}(jvar)), ')'];
+            lables_tmp_tex = [lables_tmp_tex, '_{t', num2str(M_.matched_moments{jm,2}(jvar)), '}'];
+        else
+            lables_tmp_tex = [lables_tmp_tex, '_{t}'];
+        end
+        if M_.matched_moments{jm,3}(jvar) > 1
+            lables_tmp = [lables_tmp, '^', num2str(M_.matched_moments{jm,3}(jvar))];
+            lables_tmp_tex = [lables_tmp_tex, '^{', num2str(M_.matched_moments{jm,3}(jvar)) '}'];
+        end
+        if jvar == length(M_.matched_moments{jm,1})
+            lables_tmp = [lables_tmp, ']'];
+            lables_tmp_tex = [lables_tmp_tex, ' \right]'];
+        else
+            lables_tmp = [lables_tmp, '*'];
+            lables_tmp_tex = [lables_tmp_tex, ' \times '];
+        end
+    end
+    labels{jm,1} = lables_tmp;
+    labels_TeX{jm,1} = lables_tmp_tex;
+end
+data_mat=[oo_.mom.data_moments oo_.mom.model_moments ];
 dyntable(options_mom_, title, headers, labels, data_mat, cellofchararraymaxlength(labels)+2, 10, 7);
 if options_mom_.TeX
-    lh = cellofchararraymaxlength(labels)+2;
-    labels_TeX = M_.matched_moments(:,5);
-    dyn_latex_table(M_, options_mom_, title, 'sim_corr_matrix', headers, labels_TeX, data_mat, lh, 10, 7);
+    dyn_latex_table(M_, options_mom_, title, ['comparison_moments_', options_mom_.mom.mom_method], headers, labels_TeX, data_mat, cellofchararraymaxlength(labels)+2, 10, 7);
 end
 
 if options_mom_.mode_check.status

matlab/method_of_moments/method_of_moments_objective_function.m

@@ -20,6 +20,9 @@ function [fval, info, exit_flag, df, junk1, oo_, M_, options_mom_] = method_of_m
 %   o oo_:                      structure containing the results with the following updated fields:
 %      - mom.model_moments       [numMom x 1] vector with model moments
 %      - mom.Q                   value of the quadratic form of the moment difference
+%      - mom.model_moments_params_derivs
+%                                [numMom x numParams] Jacobian matrix of derivatives of model_moments with respect to estimated parameters
+%                                (only for GMM with analytical derivatives)
 %   o M_:                       Matlab's structure describing the model
 %   o options_mom_:             structure information about all settings (specified by the user, preprocessor, and taken from global options_)
 % -------------------------------------------------------------------------
@@ -154,65 +157,55 @@ if strcmp(options_mom_.mom.mom_method,'GMM')
     end
     
     oo_.mom.model_moments = NaN(options_mom_.mom.mom_nbr,1);
-    offset = 0;
-    % First moments
-    if ~options_mom_.prefilter && isfield(options_mom_.mom.index,'E_y') && nnz(options_mom_.mom.index.E_y) > 0
-        E_y = pruned_state_space.E_y;
-        E_y_nbr = nnz(options_mom_.mom.index.E_y);
-        oo_.mom.model_moments(offset+1:E_y_nbr,1) = E_y(options_mom_.mom.index.E_y);
-        if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
-            oo_.mom.model_moments_params_derivs(offset+1:E_y_nbr,:) = pruned_state_space.dE_y(options_mom_.mom.index.E_y,:);
+    for jm = 1:size(M_.matched_moments,1)
+        % First moments
+        if ~options_mom_.prefilter && (sum(M_.matched_moments{jm,3}) == 1)
+            idx1 = (oo_.dr.obs_var == find(oo_.dr.order_var==M_.matched_moments{jm,1}) );
+            oo_.mom.model_moments(jm,1) = pruned_state_space.E_y(idx1);
+            if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
+                oo_.mom.model_moments_params_derivs(jm,:) = pruned_state_space.dE_y(idx1,:);
+            end
         end
-        offset = offset + E_y_nbr;
-    end
-    % Second moments
-    % Contemporaneous covariance
-    if isfield(options_mom_.mom.index,'E_yy') && nnz(options_mom_.mom.index.E_yy) > 0
-        if options_mom_.prefilter
-            E_yy = pruned_state_space.Var_y;
-            if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
-                dE_yy = pruned_state_space.dVar_y;
-            end            
-        else
-            E_yy = pruned_state_space.Var_y + pruned_state_space.E_y*pruned_state_space.E_y';
-            if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
-                dE_yy = pruned_state_space.dVar_y;
-                for jp=1:totparam_nbr
-                    dE_yy(:,:,jp) = dE_yy(:,:,jp) + pruned_state_space.dE_y(:,jp)*pruned_state_space.E_y' + pruned_state_space.E_y*pruned_state_space.dE_y(:,jp)';
+        % Second moments
+        if (sum(M_.matched_moments{jm,3}) == 2)
+            idx1 = (oo_.dr.obs_var == find(oo_.dr.order_var==M_.matched_moments{jm,1}(1)) );
+            idx2 = (oo_.dr.obs_var == find(oo_.dr.order_var==M_.matched_moments{jm,1}(2)) );
+            if nnz(M_.matched_moments{jm,2}) == 0
+                % Covariance
+                if options_mom_.prefilter
+                    oo_.mom.model_moments(jm,1) = pruned_state_space.Var_y(idx1,idx2);
+                    if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
+                        oo_.mom.model_moments_params_derivs(jm,:) = pruned_state_space.dVar_y(idx1,idx2,:);
+                    end
+                else
+                    oo_.mom.model_moments(jm,1) = pruned_state_space.Var_y(idx1,idx2) + pruned_state_space.E_y(idx1)*pruned_state_space.E_y(idx2)';
+                    if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
+                        for jp=1:totparam_nbr
+                            oo_.mom.model_moments_params_derivs(jm,jp) = pruned_state_space.dVar_y(idx1,idx2,jp) + pruned_state_space.dE_y(idx1,jp)*pruned_state_space.E_y(idx2)' + pruned_state_space.E_y(idx1)*pruned_state_space.dE_y(idx2,jp)';
+                        end
+                    end
+                end
+            else
+                % Autocovariance
+                lag = -M_.matched_moments{jm,2}(2); %note that leads/lags in matched_moments are transformed such that first entry is always 0 and the second is a lag
+                if options_mom_.prefilter
+                    oo_.mom.model_moments(jm,1) = pruned_state_space.Var_yi(idx1,idx2,lag);
+                    if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
+                        oo_.mom.model_moments_params_derivs(jm,:) = pruned_state_space.dVar_yi(idx1,idx2,lag,:);
+                    end
+                else
+                    oo_.mom.model_moments(jm,1) = pruned_state_space.Var_yi(idx1,idx2,lag) + pruned_state_space.E_y(idx1)*pruned_state_space.E_y(idx2)';
+                    if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
+                        for jp=1:totparam_nbr
+                            oo_.mom.model_moments_params_derivs(jm,jp) = vec( pruned_state_space.dVar_yi(idx1,idx2,lag,jp) + pruned_state_space.dE_y(idx1,jp)*pruned_state_space.E_y(idx2)' + pruned_state_space.E_y(idx1)*pruned_state_space.dE_y(idx2,jp)');
+                        end
+                    end
                 end
             end
         end
-        E_yy_nbr = nnz(tril(options_mom_.mom.index.E_yy));
-        oo_.mom.model_moments(offset+(1:E_yy_nbr),1) = E_yy(tril(options_mom_.mom.index.E_yy));
-        if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
-            oo_.mom.model_moments_params_derivs(offset+(1:E_yy_nbr),:) = reshape(dE_yy(repmat(tril(options_mom_.mom.index.E_yy),[1 1 totparam_nbr])),E_yy_nbr,totparam_nbr);
-        end
-        offset = offset + E_yy_nbr;
-    end
-    % Lead/lags covariance
-    if isfield(options_mom_.mom.index,'E_yyt') && nnz(options_mom_.mom.index.E_yyt) > 0
-        if options_mom_.prefilter
-            E_yyt = pruned_state_space.Var_yi;
-            if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
-                dE_yyt = pruned_state_space.dVar_yi;
-            end
-        else
-            E_yyt = pruned_state_space.Var_yi + repmat(pruned_state_space.E_y*pruned_state_space.E_y',[1 1 size(pruned_state_space.Var_yi,3)]);
-            if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
-                dE_yyt = pruned_state_space.dVar_yi;
-                for jp=1:totparam_nbr
-                    dE_yyt(:,:,:,jp) = dE_yyt(:,:,:,jp) + repmat(pruned_state_space.dE_y(:,jp)*pruned_state_space.E_y',[1 1 size(pruned_state_space.Var_yi,3)])...
-                                                        + repmat(pruned_state_space.E_y*pruned_state_space.dE_y(:,jp)',[1 1 size(pruned_state_space.Var_yi,3)]);
-                end
-            end
-        end
-        E_yyt_nbr = nnz(options_mom_.mom.index.E_yyt);
-        oo_.mom.model_moments(offset+(1:E_yyt_nbr),1) = E_yyt(options_mom_.mom.index.E_yyt);
-        if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
-            oo_.mom.model_moments_params_derivs(offset+(1:E_yyt_nbr),:) = reshape(dE_yyt(repmat(options_mom_.mom.index.E_yyt,[1 1 1 totparam_nbr])),E_yyt_nbr,totparam_nbr);
-        end
     end
 
+
 elseif strcmp(options_mom_.mom.mom_method,'SMM')
     %------------------------------------------------------------------------------
     % 3. Compute Moments of the model solution for normal innovations

mex/build/matlab/configure.ac

@@ -1,6 +1,6 @@
 dnl Process this file with autoconf to produce a configure script.
 
-dnl Copyright © 2009-2020 Dynare Team
+dnl Copyright © 2009-2021 Dynare Team
 dnl
 dnl This file is part of Dynare.
 dnl
@@ -18,7 +18,7 @@ dnl You should have received a copy of the GNU General Public License
 dnl along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
 AC_PREREQ([2.62])
-AC_INIT([dynare], [4.7-unstable])
+AC_INIT([dynare], [4.7.0])
 AC_CONFIG_SRCDIR([configure.ac])
 AM_INIT_AUTOMAKE([-Wall -Wno-portability foreign])
 

mex/build/octave/configure.ac

@@ -18,7 +18,7 @@ dnl You should have received a copy of the GNU General Public License
 dnl along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
 AC_PREREQ([2.62])
-AC_INIT([dynare], [4.7-unstable])
+AC_INIT([dynare], [4.7.0])
 AC_CONFIG_SRCDIR([configure.ac])
 AM_INIT_AUTOMAKE([-Wall -Wno-portability foreign])
 

preprocessor

@@ -1 +1 @@
-Subproject commit bb19d98712f2599380dfc704f98b33531d7414de
+Subproject commit 0bbba398a10abce325dd9c50b24ef3955e9c2299

tests/.gitignore

@@ -50,6 +50,7 @@ wsOct
 !/ep/mean_preserving_spread.m
 !/ep/rbcii_steady_state.m
 !/estimation/fsdat_simul.m
+!/estimation/method_of_moments/AnScho/AnScho_MoM_data_2.mat
 !/estimation/method_of_moments/RBC/RBC_MoM_steady_helper.m
 !/estimation/method_of_moments/RBC/RBC_Andreasen_Data_2.mat
 !/estimation/method_of_moments/AFVRR/AFVRR_data.mat

tests/Makefile.am

@@ -62,7 +62,13 @@ MODFILES = \
 	estimation/heteroskedastic_shocks/fs2000_het_corr.mod \
 	estimation/t_proposal/fs2000_student.mod \
 	estimation/tune_mh_jscale/fs2000.mod \
-	estimation/method_of_moments/AnScho/AnScho_MoM.mod \
+	estimation/method_of_moments/AnScho/AnScho_matched_moments.mod \
+	estimation/method_of_moments/AnScho/AnScho_GMM_estimParams0.mod \
+	estimation/method_of_moments/AnScho/AnScho_GMM_estimParams1.mod \
+	estimation/method_of_moments/AnScho/AnScho_GMM_estimParams2.mod \
+	estimation/method_of_moments/AnScho/AnScho_SMM_estimParams0.mod \
+	estimation/method_of_moments/AnScho/AnScho_SMM_estimParams1.mod \
+	estimation/method_of_moments/AnScho/AnScho_SMM_estimParams2.mod \
 	estimation/method_of_moments/RBC/RBC_MoM_Andreasen.mod \
 	estimation/method_of_moments/RBC/RBC_MoM_SMM_ME.mod \
 	estimation/method_of_moments/RBC/RBC_MoM_prefilter.mod \
@@ -71,6 +77,7 @@ MODFILES = \
 	estimation/method_of_moments/AFVRR/AFVRR_M0.mod \
 	estimation/method_of_moments/AFVRR/AFVRR_MFB.mod \
 	estimation/method_of_moments/AFVRR/AFVRR_MFB_RRA.mod \
+	estimation/system_prior_restriction/Gali_2015.mod \
 	estimation/no_init_estimation_check_first_obs/fs2000_init_check.mod \
 	estimation/example_nls.mod \
 	moments/example1_var_decomp.mod \
@@ -144,7 +151,7 @@ MODFILES = \
 	histval_initval_file/ramst_datafile.mod \
 	histval_initval_file/sim_exo_lead_lag.mod \
 	histval_initval_file/sim_exo_lead_lag_initvalf.mod \
-        ramst_normcdf_and_friends.mod \
+    ramst_normcdf_and_friends.mod \
 	ramst_vec.mod \
 	ramst_mshocks_vec.mod \
 	example1_varexo_det.mod \
@@ -288,6 +295,8 @@ MODFILES = \
 	kalman/lik_init/fs2000_lik_init_3.mod \
 	kalman/lik_init/fs2000_lik_init_4.mod \
 	kalman/lik_init/fs2000_lik_init_5.mod \
+	kalman/block/fs2000.mod \
+	kalman/block/fs2000_missing_data.mod \
 	kalman_initial_state/fs2000_smoother_only_initial_state.mod \
 	kalman_initial_state/fs2000_ns_smoother_only_initial_state.mod \
 	kalman_initial_state/fs2000_kalman_initial.mod \
@@ -1233,6 +1242,8 @@ EXTRA_DIST = \
 	estimation/fsdat_simul.m \
 	estimation/heteroskedastic_shocks/fs2000_het_model.inc \
 	estimation/heteroskedastic_shocks/fs2000_het_check.inc \
+	estimation/system_prior_restriction/Gali_2015_PC_slope.m \
+	estimation/system_prior_restriction/Gali_2015_prior_restrictions.m \
 	ep/mean_preserving_spread.m \
 	ep/rbcii_steady_state.m \
 	decision_rules/example1_results_dyn_432.mat \
@@ -1273,6 +1284,8 @@ EXTRA_DIST = \
 	lmmcp/sw-common-header.inc \
 	lmmcp/sw-common-footer.inc \
 	estimation/tune_mh_jscale/fs2000.inc \
+	estimation/method_of_moments/AnScho/AnScho_MoM_common.inc \
+	estimation/method_of_moments/AnScho/AnScho_MoM_data_2.mat \
 	estimation/method_of_moments/RBC/RBC_MoM_common.inc \
 	estimation/method_of_moments/RBC/RBC_MoM_steady_helper.m \
 	estimation/method_of_moments/AFVRR/AFVRR_common.inc \

tests/estimation/method_of_moments/AFVRR/AFVRR_M0.mod

@@ -1,303 +1,302 @@
-% DSGE model based on replication files of
-% Andreasen, Fernandez-Villaverde, Rubio-Ramirez (2018), The Pruned State-Space System for Non-Linear DSGE Models: Theory and Empirical Applications, Review of Economic Studies, 85, p. 1-49
-% Adapted for Dynare by Willi Mutschler (@wmutschl, willi@mutschler.eu), Jan 2021
-% =========================================================================
-% Copyright (C) 2021 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 <https://www.gnu.org/licenses/>.
-% =========================================================================
-
-% This is the benchmark model with no feedback M_0
-% Original code RunGMM_standardModel_RRA.m by Martin M. Andreasen, Jan 2016
-
-@#include "AFVRR_common.inc"
-
-%--------------------------------------------------------------------------
-% Parameter calibration taken from RunGMM_standardModel_RRA.m
-%--------------------------------------------------------------------------
-% fixed parameters
-INHABIT   = 1;
-PHI1      = 4;
-PHI4      = 1;
-KAPAone   = 0;
-DELTA     = 0.025;
-THETA     = 0.36;
-ETA       = 6;
-CHI       = 0;
-CONSxhr40 = 0;
-BETTAxhr  = 0;
-BETTAxhr40= 0;
-RHOD      = 0;
-GAMA      = 0.9999;
-CONSxhr20 = 0;
-
-% estimated parameters
-BETTA    = 0.999544966118000;
-B        = 0.668859504661000;
-H        = 0.342483445196000;
-PHI2     = 0.997924964981000;
-RRA      = 662.7953149595370;
-KAPAtwo  = 5.516226495551000;
-ALFA     = 0.809462321180000;
-RHOR     = 0.643873352513000;
-BETTAPAI = 1.270087844103000;
-BETTAY   = 0.031812764291000;
-MYYPS    = 1.001189151180000;
-MYZ      = 1.005286347928000;
-RHOA     = 0.743239127127000;
-RHOG     = 0.793929380230000;
-PAI      = 1.012163659169000;
-GoY      = 0.206594858866000;
-STDA     = 0.016586292524000;
-STDG     = 0.041220613851000;
-STDD     = 0.013534473123000;
-
-% endogenous parameters set via steady state, no need to initialize
-%PHIzero  = ;
-%AA       = ;
-%PHI3     = ;
-%negVf    = ;
-
-model_diagnostics;
-% Model diagnostics show that some parameters are endogenously determined
-% via the steady state, so we run steady to calibrate all parameters
-steady;
-model_diagnostics;
-% Now all parameters are determined
-
-resid;
-check;
-
-%--------------------------------------------------------------------------
-% Shock distribution
-%--------------------------------------------------------------------------
-shocks;
-var eps_a = STDA^2;
-var eps_d = STDD^2;
-var eps_g = STDG^2;
-end;
-
-%--------------------------------------------------------------------------
-% Estimated Params block - these parameters will be estimated, we
-% initialize at calibrated values
-%--------------------------------------------------------------------------
-estimated_params;
-BETTA;
-B;
-H;
-PHI2;
-RRA;
-KAPAtwo;
-ALFA;
-RHOR;
-BETTAPAI;
-BETTAY;
-MYYPS;
-MYZ;
-RHOA;
-RHOG;
-PAI;
-GoY;
-stderr eps_a;
-stderr eps_g;
-stderr eps_d;
-end;
-
-estimated_params_init(use_calibration);
-end;
-
-%--------------------------------------------------------------------------
-% Compare whether toolbox yields equivalent moments at second order
-%--------------------------------------------------------------------------
-% Note that we compare results for orderApp=1|2 and not for orderApp=3, because
-% there is a small error in the replication files of the original article in the 
-% computation of the covariance matrix of the extended innovations vector.
-% The authors have been contacted, fixed it, and report that the results 
-% change only slightly at orderApp=3 to what they report in the paper. At
-% orderApp=2 all is correct and so the following part tests whether we get 
-% the same model moments at the calibrated parameters (we do not optimize).
-% We compare it to the replication file RunGMM_standardModel_RRA.m with the
-% following settings: orderApp=1|2, seOn=0, q_lag=10, weighting=1;
-% scaled=0; optimizer=0; estimator=1; momentSet=2;
-%
-% Output of the replication files for orderApp=1
-AndreasenEtAl.Q1 = 23893.072;
-AndreasenEtAl.moments1 =[ % note that we reshuffeled to be compatible with our matched moments block
-    {[ 1]}    {'Ex'  }    {'Gr_C   '}    {'     '  }    {'0.024388'   }    {'0.023764'   }
-    {[ 2]}    {'Ex'  }    {'Gr_I   '}    {'     '  }    {'0.031046'   }    {'0.028517'   }
-    {[ 3]}    {'Ex'  }    {'Infl  ' }    {'     '  }    {'0.03757'    }    {'0.048361'   }
-    {[ 4]}    {'Ex'  }    {'r1    ' }    {'     '  }    {'0.056048'   }    {'0.073945'   }
-    {[ 5]}    {'Ex'  }    {'r40   ' }    {'     '  }    {'0.069929'   }    {'0.073945'   }
-    {[ 6]}    {'Ex'  }    {'xhr40  '}    {'     '  }    {'0.017237'   }    {'0'          }
-    {[ 7]}    {'Ex'  }    {'GoY    '}    {'     '  }    {'-1.5745'    }    {'-1.577'     }
-    {[ 8]}    {'Ex'  }    {'hours  '}    {'     '  }    {'-0.043353'  }    {'-0.042861'  }
-    {[ 9]}    {'Exx' }    {'Gr_C   '}    {'Gr_C   '}    {'0.0013159'  }    {'0.0011816'  }
-    {[17]}    {'Exx' }    {'Gr_C   '}    {'Gr_I   '}    {'0.0021789'  }    {'0.0016052'  }
-    {[18]}    {'Exx' }    {'Gr_C   '}    {'Infl  ' }    {'0.00067495' }    {'0.00090947' }
-    {[19]}    {'Exx' }    {'Gr_C   '}    {'r1    ' }    {'0.0011655'  }    {'0.0016016'  }
-    {[20]}    {'Exx' }    {'Gr_C   '}    {'r40   ' }    {'0.0015906'  }    {'0.0017076'  }
-    {[21]}    {'Exx' }    {'Gr_C   '}    {'xhr40  '}    {'0.0020911'  }    {'0.0013997'  }
-    {[10]}    {'Exx' }    {'Gr_I   '}    {'Gr_I   '}    {'0.0089104'  }    {'0.0055317'  }
-    {[22]}    {'Exx' }    {'Gr_I   '}    {'Infl  ' }    {'0.00063139' }    {'0.00050106' }
-    {[23]}    {'Exx' }    {'Gr_I   '}    {'r1    ' }    {'0.0011031'  }    {'0.0018178'  }
-    {[24]}    {'Exx' }    {'Gr_I   '}    {'r40   ' }    {'0.0018445'  }    {'0.0020186'  }
-    {[25]}    {'Exx' }    {'Gr_I   '}    {'xhr40  '}    {'0.00095556' }    {'0.0064471'  }
-    {[11]}    {'Exx' }    {'Infl  ' }    {'Infl  ' }    {'0.0020268'  }    {'0.0030519'  }
-    {[26]}    {'Exx' }    {'Infl  ' }    {'r1    ' }    {'0.0025263'  }    {'0.0042181'  }
-    {[27]}    {'Exx' }    {'Infl  ' }    {'r40   ' }    {'0.0029126'  }    {'0.0039217'  }
-    {[28]}    {'Exx' }    {'Infl  ' }    {'xhr40  '}    {'-0.00077101'}    {'-0.0019975' }
-    {[12]}    {'Exx' }    {'r1    ' }    {'r1    ' }    {'0.0038708'  }    {'0.0061403'  }
-    {[29]}    {'Exx' }    {'r1    ' }    {'r40   ' }    {'0.0044773'  }    {'0.0058343'  }
-    {[30]}    {'Exx' }    {'r1    ' }    {'xhr40  '}    {'-0.00048202'}    {'-0.00089501'}
-    {[13]}    {'Exx' }    {'r40   ' }    {'r40   ' }    {'0.0054664'  }    {'0.0056883'  }
-    {[31]}    {'Exx' }    {'r40   ' }    {'xhr40  '}    {'0.00053864' }    {'-0.00041184'}
-    {[14]}    {'Exx' }    {'xhr40  '}    {'xhr40  '}    {'0.053097'   }    {'0.016255'   }
-    {[15]}    {'Exx' }    {'GoY    '}    {'GoY    '}    {'2.4863'     }    {'2.4919'     }
-    {[16]}    {'Exx' }    {'hours  '}    {'hours  '}    {'0.0018799'  }    {'0.0018384'  }
-    {[32]}    {'Exx1'}    {'Gr_C   '}    {'Gr_C   '}    {'0.00077917' }    {'0.00065543' }
-    {[33]}    {'Exx1'}    {'Gr_I   '}    {'Gr_I   '}    {'0.0050104'  }    {'0.0033626'  }
-    {[34]}    {'Exx1'}    {'Infl  ' }    {'Infl  ' }    {'0.0019503'  }    {'0.0029033'  }
-    {[35]}    {'Exx1'}    {'r1    ' }    {'r1    ' }    {'0.0038509'  }    {'0.006112'   }
-    {[36]}    {'Exx1'}    {'r40   ' }    {'r40   ' }    {'0.0054699'  }    {'0.005683'   }
-    {[37]}    {'Exx1'}    {'xhr40  '}    {'xhr40  '}    {'-0.00098295'}    {'3.3307e-16' }
-    {[38]}    {'Exx1'}    {'GoY    '}    {'GoY    '}    {'2.4868'     }    {'2.4912'     }
-    {[39]}    {'Exx1'}    {'hours  '}    {'hours  '}    {'0.0018799'  }    {'0.0018378'  }
-];
-
-% Output of the replication files for orderApp=2
-AndreasenEtAl.Q2 = 65.8269;
-AndreasenEtAl.moments2 =[ % note that we reshuffeled to be compatible with our matched moments block
-    {[ 1]}    {'Ex'  }    {'Gr_C   '}    {'     '  }    {'0.024388'   }    {'0.023764'  }
-    {[ 2]}    {'Ex'  }    {'Gr_I   '}    {'     '  }    {'0.031046'   }    {'0.028517'  }
-    {[ 3]}    {'Ex'  }    {'Infl  ' }    {'     '  }    {'0.03757'    }    {'0.034882'  }
-    {[ 4]}    {'Ex'  }    {'r1    ' }    {'     '  }    {'0.056048'   }    {'0.056542'  }
-    {[ 5]}    {'Ex'  }    {'r40   ' }    {'     '  }    {'0.069929'   }    {'0.070145'  }
-    {[ 6]}    {'Ex'  }    {'xhr40  '}    {'     '  }    {'0.017237'   }    {'0.020825'  }
-    {[ 7]}    {'Ex'  }    {'GoY    '}    {'     '  }    {'-1.5745'    }    {'-1.5748'   }
-    {[ 8]}    {'Ex'  }    {'hours  '}    {'     '  }    {'-0.043353'  }    {'-0.04335'  }
-    {[ 9]}    {'Exx' }    {'Gr_C   '}    {'Gr_C   '}    {'0.0013159'  }    {'0.001205'  }
-    {[17]}    {'Exx' }    {'Gr_C   '}    {'Gr_I   '}    {'0.0021789'  }    {'0.0016067' }
-    {[18]}    {'Exx' }    {'Gr_C   '}    {'Infl  ' }    {'0.00067495' }    {'0.00059406'}
-    {[19]}    {'Exx' }    {'Gr_C   '}    {'r1    ' }    {'0.0011655'  }    {'0.0011949' }
-    {[20]}    {'Exx' }    {'Gr_C   '}    {'r40   ' }    {'0.0015906'  }    {'0.0016104' }
-    {[21]}    {'Exx' }    {'Gr_C   '}    {'xhr40  '}    {'0.0020911'  }    {'0.0020245' }
-    {[10]}    {'Exx' }    {'Gr_I   '}    {'Gr_I   '}    {'0.0089104'  }    {'0.0060254' }
-    {[22]}    {'Exx' }    {'Gr_I   '}    {'Infl  ' }    {'0.00063139' }    {'8.3563e-05'}
-    {[23]}    {'Exx' }    {'Gr_I   '}    {'r1    ' }    {'0.0011031'  }    {'0.0013176' }
-    {[24]}    {'Exx' }    {'Gr_I   '}    {'r40   ' }    {'0.0018445'  }    {'0.0019042' }
-    {[25]}    {'Exx' }    {'Gr_I   '}    {'xhr40  '}    {'0.00095556' }    {'0.0064261' }
-    {[11]}    {'Exx' }    {'Infl  ' }    {'Infl  ' }    {'0.0020268'  }    {'0.0020735' }
-    {[26]}    {'Exx' }    {'Infl  ' }    {'r1    ' }    {'0.0025263'  }    {'0.0027621' }
-    {[27]}    {'Exx' }    {'Infl  ' }    {'r40   ' }    {'0.0029126'  }    {'0.0029257' }
-    {[28]}    {'Exx' }    {'Infl  ' }    {'xhr40  '}    {'-0.00077101'}    {'-0.0012165'}
-    {[12]}    {'Exx' }    {'r1    ' }    {'r1    ' }    {'0.0038708'  }    {'0.0040235' }
-    {[29]}    {'Exx' }    {'r1    ' }    {'r40   ' }    {'0.0044773'  }    {'0.0044702' }
-    {[30]}    {'Exx' }    {'r1    ' }    {'xhr40  '}    {'-0.00048202'}    {'0.00030542'}
-    {[13]}    {'Exx' }    {'r40   ' }    {'r40   ' }    {'0.0054664'  }    {'0.0052718' }
-    {[31]}    {'Exx' }    {'r40   ' }    {'xhr40  '}    {'0.00053864' }    {'0.0010045' }
-    {[14]}    {'Exx' }    {'xhr40  '}    {'xhr40  '}    {'0.053097'   }    {'0.018416'  }
-    {[15]}    {'Exx' }    {'GoY    '}    {'GoY    '}    {'2.4863'     }    {'2.4853'    }
-    {[16]}    {'Exx' }    {'hours  '}    {'hours  '}    {'0.0018799'  }    {'0.0018806' }
-    {[32]}    {'Exx1'}    {'Gr_C   '}    {'Gr_C   '}    {'0.00077917' }    {'0.00067309'}
-    {[33]}    {'Exx1'}    {'Gr_I   '}    {'Gr_I   '}    {'0.0050104'  }    {'0.0033293' }
-    {[34]}    {'Exx1'}    {'Infl  ' }    {'Infl  ' }    {'0.0019503'  }    {'0.0019223' }
-    {[35]}    {'Exx1'}    {'r1    ' }    {'r1    ' }    {'0.0038509'  }    {'0.0039949' }
-    {[36]}    {'Exx1'}    {'r40   ' }    {'r40   ' }    {'0.0054699'  }    {'0.0052659' }
-    {[37]}    {'Exx1'}    {'xhr40  '}    {'xhr40  '}    {'-0.00098295'}    {'0.0004337' }
-    {[38]}    {'Exx1'}    {'GoY    '}    {'GoY    '}    {'2.4868'     }    {'2.4846'    }
-    {[39]}    {'Exx1'}    {'hours  '}    {'hours  '}    {'0.0018799'  }    {'0.00188'   }
-];
-
-@#for orderApp in 1:2
-
-method_of_moments(
-          mom_method = GMM                    % method of moments method; possible values: GMM|SMM
-        , datafile   = 'AFVRR_data.mat'       % name of filename with data
-        , bartlett_kernel_lag = 10          % bandwith in optimal weighting matrix
-        , order = @{orderApp}                 % order of Taylor approximation in perturbation
-        , pruning                             % use pruned state space system at higher-order
-        % , verbose                           % display and store intermediate estimation results
-        , weighting_matrix = ['DIAGONAL']      % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename        
-        % , TeX                               % print TeX tables and graphics
-        % Optimization options that can be set by the user in the mod file, otherwise default values are provided        
-        %, huge_number=1D10                   % value for replacing the infinite bounds on parameters by finite numbers. Used by some optimizers for numerical reasons
-        , mode_compute = 0                    % specifies the optimizer for minimization of moments distance, note that by default there is a new optimizer        
-        , optim = ('TolFun', 1e-6
-                   ,'TolX', 1e-6
-                   ,'MaxIter', 3000
-                   ,'MaxFunEvals', 1D6
-                   ,'UseParallel' , 1
-                   %,'Jacobian' , 'on'
-                  )    % a list of NAME and VALUE pairs to set options for the optimization routines. Available options depend on mode_compute
-        %, silent_optimizer                  % run minimization of moments distance silently without displaying results or saving files in between
-        %, analytic_standard_errors
-        , se_tolx=1e-10
-);
-
-% Check results
-
-fprintf('****************************************************************\n')
-fprintf('Compare Results for perturbation order @{orderApp}\n')
-fprintf('****************************************************************\n')
-dev_Q            = AndreasenEtAl.Q@{orderApp} - oo_.mom.Q;
-dev_datamoments  = str2double(AndreasenEtAl.moments@{orderApp}(:,5)) - oo_.mom.data_moments;
-dev_modelmoments = str2double(AndreasenEtAl.moments@{orderApp}(:,6)) - oo_.mom.model_moments;
-
-% There is no table command in Octave
-% The table command also crashes on MATLAB R2014a because it does not like variable names
-if ~isoctave && ~matlab_ver_less_than('8.4')
-table([AndreasenEtAl.Q@{orderApp} ; str2double(AndreasenEtAl.moments@{orderApp}(:,5)) ; str2double(AndreasenEtAl.moments@{orderApp}(:,6))],...
-      [oo_.mom.Q                  ; oo_.mom.data_moments                              ; oo_.mom.model_moments                            ],...
-      [dev_Q                      ; dev_datamoments                                   ; dev_modelmoments                                 ],...
-      'VariableNames', {'Andreasen et al', 'Dynare', 'dev'},...
-      'RowNames', ['Q'; strcat('Data_', M_.matched_moments(:,4)); strcat('Model_', M_.matched_moments(:,4))])
-end
-
-if norm(dev_modelmoments)> 1e-4
-    error('Something wrong in the computation of moments at order @{orderApp}')
-end
-
-@#endfor
-
-%--------------------------------------------------------------------------
-% Replicate estimation at orderApp=3
-%--------------------------------------------------------------------------
-@#ifdef DoEstimation
-method_of_moments(
-          mom_method = GMM                    % method of moments method; possible values: GMM|SMM
-        , datafile   = 'AFVRR_data.mat'       % name of filename with data
-        , bartlett_kernel_lag = 10            % bandwith in optimal weighting matrix
-        , order = 3                           % order of Taylor approximation in perturbation
-        , pruning                             % use pruned state space system at higher-order
-        % , verbose                           % display and store intermediate estimation results
-        , weighting_matrix = ['DIAGONAL', 'OPTIMAL']      % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename        
-        % , TeX                               % print TeX tables and graphics
-        % Optimization options that can be set by the user in the mod file, otherwise default values are provided        
-        %, huge_number=1D10                   % value for replacing the infinite bounds on parameters by finite numbers. Used by some optimizers for numerical reasons
-        , mode_compute = 13                   % specifies the optimizer for minimization of moments distance, note that by default there is a new optimizer
-        , additional_optimizer_steps = [13]
-        , optim = ('TolFun', 1e-6
-                   ,'TolX', 1e-6
-                   ,'MaxIter', 3000
-                   ,'MaxFunEvals', 1D6
-                   ,'UseParallel' , 1
-                   %,'Jacobian' , 'on'
-                  )    % a list of NAME and VALUE pairs to set options for the optimization routines. Available options depend on mode_compute
-        %, silent_optimizer                  % run minimization of moments distance silently without displaying results or saving files in between
-        %, analytic_standard_errors
-        , se_tolx=1e-10
-);
-@#endif
+% DSGE model based on replication files of
+% Andreasen, Fernandez-Villaverde, Rubio-Ramirez (2018), The Pruned State-Space System for Non-Linear DSGE Models: Theory and Empirical Applications, Review of Economic Studies, 85, p. 1-49
+% Adapted for Dynare by Willi Mutschler (@wmutschl, willi@mutschler.eu), Jan 2021
+% =========================================================================
+% Copyright (C) 2021 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 <https://www.gnu.org/licenses/>.
+% =========================================================================
+
+% This is the benchmark model with no feedback M_0
+% Original code RunGMM_standardModel_RRA.m by Martin M. Andreasen, Jan 2016
+
+@#include "AFVRR_common.inc"
+
+%--------------------------------------------------------------------------
+% Parameter calibration taken from RunGMM_standardModel_RRA.m
+%--------------------------------------------------------------------------
+% fixed parameters
+INHABIT   = 1;
+PHI1      = 4;
+PHI4      = 1;
+KAPAone   = 0;
+DELTA     = 0.025;
+THETA     = 0.36;
+ETA       = 6;
+CHI       = 0;
+CONSxhr40 = 0;
+BETTAxhr  = 0;
+BETTAxhr40= 0;
+RHOD      = 0;
+GAMA      = 0.9999;
+CONSxhr20 = 0;
+
+% estimated parameters
+BETTA    = 0.999544966118000;
+B        = 0.668859504661000;
+H        = 0.342483445196000;
+PHI2     = 0.997924964981000;
+RRA      = 662.7953149595370;
+KAPAtwo  = 5.516226495551000;
+ALFA     = 0.809462321180000;
+RHOR     = 0.643873352513000;
+BETTAPAI = 1.270087844103000;
+BETTAY   = 0.031812764291000;
+MYYPS    = 1.001189151180000;
+MYZ      = 1.005286347928000;
+RHOA     = 0.743239127127000;
+RHOG     = 0.793929380230000;
+PAI      = 1.012163659169000;
+GoY      = 0.206594858866000;
+STDA     = 0.016586292524000;
+STDG     = 0.041220613851000;
+STDD     = 0.013534473123000;
+
+% endogenous parameters set via steady state, no need to initialize
+%PHIzero  = ;
+%AA       = ;
+%PHI3     = ;
+%negVf    = ;
+
+model_diagnostics;
+% Model diagnostics show that some parameters are endogenously determined
+% via the steady state, so we run steady to calibrate all parameters
+steady;
+model_diagnostics;
+% Now all parameters are determined
+
+resid;
+check;
+
+%--------------------------------------------------------------------------
+% Shock distribution
+%--------------------------------------------------------------------------
+shocks;
+var eps_a = STDA^2;
+var eps_d = STDD^2;
+var eps_g = STDG^2;
+end;
+
+%--------------------------------------------------------------------------
+% Estimated Params block - these parameters will be estimated, we
+% initialize at calibrated values
+%--------------------------------------------------------------------------
+estimated_params;
+BETTA;
+B;
+H;
+PHI2;
+RRA;
+KAPAtwo;
+ALFA;
+RHOR;
+BETTAPAI;
+BETTAY;
+MYYPS;
+MYZ;
+RHOA;
+RHOG;
+PAI;
+GoY;
+stderr eps_a;
+stderr eps_g;
+stderr eps_d;
+end;
+
+estimated_params_init(use_calibration);
+end;
+
+%--------------------------------------------------------------------------
+% Compare whether toolbox yields equivalent moments at second order
+%--------------------------------------------------------------------------
+% Note that we compare results for orderApp=1|2 and not for orderApp=3, because
+% there is a small error in the replication files of the original article in the 
+% computation of the covariance matrix of the extended innovations vector.
+% The authors have been contacted, fixed it, and report that the results 
+% change only slightly at orderApp=3 to what they report in the paper. At
+% orderApp=2 all is correct and so the following part tests whether we get 
+% the same model moments at the calibrated parameters (we do not optimize).
+% We compare it to the replication file RunGMM_standardModel_RRA.m with the
+% following settings: orderApp=1|2, seOn=0, q_lag=10, weighting=1;
+% scaled=0; optimizer=0; estimator=1; momentSet=2;
+%
+% Output of the replication files for orderApp=1
+AndreasenEtAl.Q1 = 23893.072;
+AndreasenEtAl.moments1 =[ % note that we reshuffeled to be compatible with our matched moments block
+    {[ 1]}    {'Ex'  }    {'Gr_C   '}    {'     '  }    {'0.024388'   }    {'0.023764'   }
+    {[ 2]}    {'Ex'  }    {'Gr_I   '}    {'     '  }    {'0.031046'   }    {'0.028517'   }
+    {[ 3]}    {'Ex'  }    {'Infl  ' }    {'     '  }    {'0.03757'    }    {'0.048361'   }
+    {[ 4]}    {'Ex'  }    {'r1    ' }    {'     '  }    {'0.056048'   }    {'0.073945'   }
+    {[ 5]}    {'Ex'  }    {'r40   ' }    {'     '  }    {'0.069929'   }    {'0.073945'   }
+    {[ 6]}    {'Ex'  }    {'xhr40  '}    {'     '  }    {'0.017237'   }    {'0'          }
+    {[ 7]}    {'Ex'  }    {'GoY    '}    {'     '  }    {'-1.5745'    }    {'-1.577'     }
+    {[ 8]}    {'Ex'  }    {'hours  '}    {'     '  }    {'-0.043353'  }    {'-0.042861'  }
+    {[ 9]}    {'Exx' }    {'Gr_C   '}    {'Gr_C   '}    {'0.0013159'  }    {'0.0011816'  }
+    {[10]}    {'Exx' }    {'Gr_C   '}    {'Gr_I   '}    {'0.0021789'  }    {'0.0016052'  }
+    {[11]}    {'Exx' }    {'Gr_C   '}    {'Infl  ' }    {'0.00067495' }    {'0.00090947' }
+    {[12]}    {'Exx' }    {'Gr_C   '}    {'r1    ' }    {'0.0011655'  }    {'0.0016016'  }
+    {[13]}    {'Exx' }    {'Gr_C   '}    {'r40   ' }    {'0.0015906'  }    {'0.0017076'  }
+    {[14]}    {'Exx' }    {'Gr_C   '}    {'xhr40  '}    {'0.0020911'  }    {'0.0013997'  }
+    {[15]}    {'Exx' }    {'Gr_I   '}    {'Gr_I   '}    {'0.0089104'  }    {'0.0055317'  }
+    {[16]}    {'Exx' }    {'Gr_I   '}    {'Infl  ' }    {'0.00063139' }    {'0.00050106' }
+    {[17]}    {'Exx' }    {'Gr_I   '}    {'r1    ' }    {'0.0011031'  }    {'0.0018178'  }
+    {[18]}    {'Exx' }    {'Gr_I   '}    {'r40   ' }    {'0.0018445'  }    {'0.0020186'  }
+    {[19]}    {'Exx' }    {'Gr_I   '}    {'xhr40  '}    {'0.00095556' }    {'0.0064471'  }
+    {[20]}    {'Exx' }    {'Infl  ' }    {'Infl  ' }    {'0.0020268'  }    {'0.0030519'  }
+    {[21]}    {'Exx' }    {'Infl  ' }    {'r1    ' }    {'0.0025263'  }    {'0.0042181'  }
+    {[22]}    {'Exx' }    {'Infl  ' }    {'r40   ' }    {'0.0029126'  }    {'0.0039217'  }
+    {[23]}    {'Exx' }    {'Infl  ' }    {'xhr40  '}    {'-0.00077101'}    {'-0.0019975' }
+    {[24]}    {'Exx' }    {'r1    ' }    {'r1    ' }    {'0.0038708'  }    {'0.0061403'  }
+    {[25]}    {'Exx' }    {'r1    ' }    {'r40   ' }    {'0.0044773'  }    {'0.0058343'  }
+    {[26]}    {'Exx' }    {'r1    ' }    {'xhr40  '}    {'-0.00048202'}    {'-0.00089501'}
+    {[27]}    {'Exx' }    {'r40   ' }    {'r40   ' }    {'0.0054664'  }    {'0.0056883'  }
+    {[28]}    {'Exx' }    {'r40   ' }    {'xhr40  '}    {'0.00053864' }    {'-0.00041184'}
+    {[29]}    {'Exx' }    {'xhr40  '}    {'xhr40  '}    {'0.053097'   }    {'0.016255'   }
+    {[30]}    {'Exx' }    {'GoY    '}    {'GoY    '}    {'2.4863'     }    {'2.4919'     }
+    {[31]}    {'Exx' }    {'hours  '}    {'hours  '}    {'0.0018799'  }    {'0.0018384'  }
+    {[32]}    {'Exx1'}    {'Gr_C   '}    {'Gr_C   '}    {'0.00077917' }    {'0.00065543' }
+    {[33]}    {'Exx1'}    {'Gr_I   '}    {'Gr_I   '}    {'0.0050104'  }    {'0.0033626'  }
+    {[34]}    {'Exx1'}    {'Infl  ' }    {'Infl  ' }    {'0.0019503'  }    {'0.0029033'  }
+    {[35]}    {'Exx1'}    {'r1    ' }    {'r1    ' }    {'0.0038509'  }    {'0.006112'   }
+    {[36]}    {'Exx1'}    {'r40   ' }    {'r40   ' }    {'0.0054699'  }    {'0.005683'   }
+    {[37]}    {'Exx1'}    {'xhr40  '}    {'xhr40  '}    {'-0.00098295'}    {'3.3307e-16' }
+    {[38]}    {'Exx1'}    {'GoY    '}    {'GoY    '}    {'2.4868'     }    {'2.4912'     }
+    {[39]}    {'Exx1'}    {'hours  '}    {'hours  '}    {'0.0018799'  }    {'0.0018378'  }
+];
+
+% Output of the replication files for orderApp=2
+AndreasenEtAl.Q2 = 65.8269;
+AndreasenEtAl.moments2 =[ % note that we reshuffeled to be compatible with our matched moments block
+    {[ 1]}    {'Ex'  }    {'Gr_C   '}    {'     '  }    {'0.024388'   }    {'0.023764'  }
+    {[ 2]}    {'Ex'  }    {'Gr_I   '}    {'     '  }    {'0.031046'   }    {'0.028517'  }
+    {[ 3]}    {'Ex'  }    {'Infl  ' }    {'     '  }    {'0.03757'    }    {'0.034882'  }
+    {[ 4]}    {'Ex'  }    {'r1    ' }    {'     '  }    {'0.056048'   }    {'0.056542'  }
+    {[ 5]}    {'Ex'  }    {'r40   ' }    {'     '  }    {'0.069929'   }    {'0.070145'  }
+    {[ 6]}    {'Ex'  }    {'xhr40  '}    {'     '  }    {'0.017237'   }    {'0.020825'  }
+    {[ 7]}    {'Ex'  }    {'GoY    '}    {'     '  }    {'-1.5745'    }    {'-1.5748'   }
+    {[ 8]}    {'Ex'  }    {'hours  '}    {'     '  }    {'-0.043353'  }    {'-0.04335'  }
+    {[ 9]}    {'Exx' }    {'Gr_C   '}    {'Gr_C   '}    {'0.0013159'  }    {'0.001205'  }
+    {[10]}    {'Exx' }    {'Gr_C   '}    {'Gr_I   '}    {'0.0021789'  }    {'0.0016067' }
+    {[11]}    {'Exx' }    {'Gr_C   '}    {'Infl  ' }    {'0.00067495' }    {'0.00059406'}
+    {[12]}    {'Exx' }    {'Gr_C   '}    {'r1    ' }    {'0.0011655'  }    {'0.0011949' }
+    {[13]}    {'Exx' }    {'Gr_C   '}    {'r40   ' }    {'0.0015906'  }    {'0.0016104' }
+    {[14]}    {'Exx' }    {'Gr_C   '}    {'xhr40  '}    {'0.0020911'  }    {'0.0020245' }
+    {[15]}    {'Exx' }    {'Gr_I   '}    {'Gr_I   '}    {'0.0089104'  }    {'0.0060254' }
+    {[16]}    {'Exx' }    {'Gr_I   '}    {'Infl  ' }    {'0.00063139' }    {'8.3563e-05'}
+    {[17]}    {'Exx' }    {'Gr_I   '}    {'r1    ' }    {'0.0011031'  }    {'0.0013176' }
+    {[18]}    {'Exx' }    {'Gr_I   '}    {'r40   ' }    {'0.0018445'  }    {'0.0019042' }
+    {[19]}    {'Exx' }    {'Gr_I   '}    {'xhr40  '}    {'0.00095556' }    {'0.0064261' }
+    {[20]}    {'Exx' }    {'Infl  ' }    {'Infl  ' }    {'0.0020268'  }    {'0.0020735' }
+    {[21]}    {'Exx' }    {'Infl  ' }    {'r1    ' }    {'0.0025263'  }    {'0.0027621' }
+    {[22]}    {'Exx' }    {'Infl  ' }    {'r40   ' }    {'0.0029126'  }    {'0.0029257' }
+    {[23]}    {'Exx' }    {'Infl  ' }    {'xhr40  '}    {'-0.00077101'}    {'-0.0012165'}
+    {[24]}    {'Exx' }    {'r1    ' }    {'r1    ' }    {'0.0038708'  }    {'0.0040235' }
+    {[25]}    {'Exx' }    {'r1    ' }    {'r40   ' }    {'0.0044773'  }    {'0.0044702' }
+    {[26]}    {'Exx' }    {'r1    ' }    {'xhr40  '}    {'-0.00048202'}    {'0.00030542'}
+    {[27]}    {'Exx' }    {'r40   ' }    {'r40   ' }    {'0.0054664'  }    {'0.0052718' }
+    {[28]}    {'Exx' }    {'r40   ' }    {'xhr40  '}    {'0.00053864' }    {'0.0010045' }
+    {[29]}    {'Exx' }    {'xhr40  '}    {'xhr40  '}    {'0.053097'   }    {'0.018416'  }
+    {[30]}    {'Exx' }    {'GoY    '}    {'GoY    '}    {'2.4863'     }    {'2.4853'    }
+    {[31]}    {'Exx' }    {'hours  '}    {'hours  '}    {'0.0018799'  }    {'0.0018806' }
+    {[32]}    {'Exx1'}    {'Gr_C   '}    {'Gr_C   '}    {'0.00077917' }    {'0.00067309'}
+    {[33]}    {'Exx1'}    {'Gr_I   '}    {'Gr_I   '}    {'0.0050104'  }    {'0.0033293' }
+    {[34]}    {'Exx1'}    {'Infl  ' }    {'Infl  ' }    {'0.0019503'  }    {'0.0019223' }
+    {[35]}    {'Exx1'}    {'r1    ' }    {'r1    ' }    {'0.0038509'  }    {'0.0039949' }
+    {[36]}    {'Exx1'}    {'r40   ' }    {'r40   ' }    {'0.0054699'  }    {'0.0052659' }
+    {[37]}    {'Exx1'}    {'xhr40  '}    {'xhr40  '}    {'-0.00098295'}    {'0.0004337' }
+    {[38]}    {'Exx1'}    {'GoY    '}    {'GoY    '}    {'2.4868'     }    {'2.4846'    }
+    {[39]}    {'Exx1'}    {'hours  '}    {'hours  '}    {'0.0018799'  }    {'0.00188'   }
+];
+
+@#for orderApp in 1:2
+
+method_of_moments(
+          mom_method = GMM                    % method of moments method; possible values: GMM|SMM
+        , datafile   = 'AFVRR_data.mat'       % name of filename with data
+        , bartlett_kernel_lag = 10          % bandwith in optimal weighting matrix
+        , order = @{orderApp}                 % order of Taylor approximation in perturbation
+        , pruning                             % use pruned state space system at higher-order
+        % , verbose                           % display and store intermediate estimation results
+        , weighting_matrix = ['DIAGONAL']      % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename        
+        % , TeX                               % print TeX tables and graphics
+        % Optimization options that can be set by the user in the mod file, otherwise default values are provided        
+        %, huge_number=1D10                   % value for replacing the infinite bounds on parameters by finite numbers. Used by some optimizers for numerical reasons
+        , mode_compute = 0                    % specifies the optimizer for minimization of moments distance, note that by default there is a new optimizer        
+        , optim = ('TolFun', 1e-6
+                   ,'TolX', 1e-6
+                   ,'MaxIter', 3000
+                   ,'MaxFunEvals', 1D6
+                   ,'UseParallel' , 1
+                   %,'Jacobian' , 'on'
+                  )    % a list of NAME and VALUE pairs to set options for the optimization routines. Available options depend on mode_compute
+        %, silent_optimizer                  % run minimization of moments distance silently without displaying results or saving files in between
+        %, analytic_standard_errors
+        , se_tolx=1e-10
+);
+
+% Check results
+
+fprintf('****************************************************************\n')
+fprintf('Compare Results for perturbation order @{orderApp}\n')
+fprintf('****************************************************************\n')
+dev_Q            = AndreasenEtAl.Q@{orderApp} - oo_.mom.Q;
+dev_datamoments  = str2double(AndreasenEtAl.moments@{orderApp}(:,5)) - oo_.mom.data_moments;
+dev_modelmoments = str2double(AndreasenEtAl.moments@{orderApp}(:,6)) - oo_.mom.model_moments;
+
+% There is no table command in Octave
+% The table command also crashes on MATLAB R2014a because it does not like variable names
+if ~isoctave && ~matlab_ver_less_than('8.4')
+table([AndreasenEtAl.Q@{orderApp} ; str2double(AndreasenEtAl.moments@{orderApp}(:,5)) ; str2double(AndreasenEtAl.moments@{orderApp}(:,6))],...
+      [oo_.mom.Q                  ; oo_.mom.data_moments                              ; oo_.mom.model_moments                            ],...
+      [dev_Q                      ; dev_datamoments                                   ; dev_modelmoments                                 ],...
+      'VariableNames', {'Andreasen et al', 'Dynare', 'dev'})
+end
+
+if norm(dev_modelmoments)> 1e-4
+    error('Something wrong in the computation of moments at order @{orderApp}')
+end
+
+@#endfor
+
+%--------------------------------------------------------------------------
+% Replicate estimation at orderApp=3
+%--------------------------------------------------------------------------
+@#ifdef DoEstimation
+method_of_moments(
+          mom_method = GMM                    % method of moments method; possible values: GMM|SMM
+        , datafile   = 'AFVRR_data.mat'       % name of filename with data
+        , bartlett_kernel_lag = 10            % bandwith in optimal weighting matrix
+        , order = 3                           % order of Taylor approximation in perturbation
+        , pruning                             % use pruned state space system at higher-order
+        % , verbose                           % display and store intermediate estimation results
+        , weighting_matrix = ['DIAGONAL', 'OPTIMAL']      % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename        
+        % , TeX                               % print TeX tables and graphics
+        % Optimization options that can be set by the user in the mod file, otherwise default values are provided        
+        %, huge_number=1D10                   % value for replacing the infinite bounds on parameters by finite numbers. Used by some optimizers for numerical reasons
+        , mode_compute = 13                   % specifies the optimizer for minimization of moments distance, note that by default there is a new optimizer
+        , additional_optimizer_steps = [13]
+        , optim = ('TolFun', 1e-6
+                   ,'TolX', 1e-6
+                   ,'MaxIter', 3000
+                   ,'MaxFunEvals', 1D6
+                   ,'UseParallel' , 1
+                   %,'Jacobian' , 'on'
+                  )    % a list of NAME and VALUE pairs to set options for the optimization routines. Available options depend on mode_compute
+        %, silent_optimizer                  % run minimization of moments distance silently without displaying results or saving files in between
+        %, analytic_standard_errors
+        , se_tolx=1e-10
+);
+@#endif

tests/estimation/method_of_moments/AFVRR/AFVRR_MFB.mod

@@ -1,304 +1,303 @@
-% DSGE model based on replication files of
-% Andreasen, Fernandez-Villaverde, Rubio-Ramirez (2018), The Pruned State-Space System for Non-Linear DSGE Models: Theory and Empirical Applications, Review of Economic Studies, 85, p. 1-49
-% Adapted for Dynare by Willi Mutschler (@wmutschl, willi@mutschler.eu), Jan 2021
-% =========================================================================
-% Copyright (C) 2021 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 <https://www.gnu.org/licenses/>.
-% =========================================================================
-
-% This is the model with Feedback M_FB
-% Original code RunGMM_Feedback_estim_RRA.m by Martin M. Andreasen, Jan 2016
-
-@#include "AFVRR_common.inc"
-
-%--------------------------------------------------------------------------
-% Parameter calibration taken from RunGMM_Feedback_estim_RRA.m
-%--------------------------------------------------------------------------
-% fixed parameters
-INHABIT   = 1;
-PHI1      = 4;
-PHI4      = 1;
-KAPAone   = 0;
-DELTA     = 0.025;
-THETA     = 0.36;
-ETA       = 6;
-CHI       = 0;
-BETTAxhr  = 0;
-BETTAxhr40= 0;
-RHOD      = 0;
-GAMA      = 0.9999;
-CONSxhr20 = 0;
-
-% estimated parameters
-BETTA     = 0.997007023687000;
-B         = 0.692501768577000;
-H         = 0.339214495653000;
-PHI2      = 0.688555040951000;
-RRA       = 24.346514272871001;
-KAPAtwo   = 10.018421876923000;
-ALFA      = 0.792507553312000;
-RHOR      = 0.849194030384000;
-BETTAPAI  = 2.060579322980000;
-BETTAY    = 0.220573712342000;
-MYYPS     = 1.001016690133000;
-MYZ       = 1.005356313981000;
-RHOA      = 0.784141391843000;
-RHOG      = 0.816924540497000;
-PAI       = 1.011924196487000;
-CONSxhr40 = 0.878774662208000;
-GoY       = 0.207110300602000;
-STDA      = 0.013024450606000;
-STDG      = 0.051049871928000;
-STDD      = 0.008877423780000;
-
-% endogenous parameters set via steady state, no need to initialize
-%PHIzero  = ;
-%AA       = ;
-%PHI3     = ;
-%negVf    = ;
-
-model_diagnostics;
-% Model diagnostics show that some parameters are endogenously determined
-% via the steady state, so we run steady to calibrate all parameters
-steady;
-model_diagnostics;
-% Now all parameters are determined
-
-resid;
-check;
-
-%--------------------------------------------------------------------------
-% Shock distribution
-%--------------------------------------------------------------------------
-shocks;
-var eps_a = STDA^2;
-var eps_d = STDD^2;
-var eps_g = STDG^2;
-end;
-
-%--------------------------------------------------------------------------
-% Estimated Params block - these parameters will be estimated, we
-% initialize at calibrated values
-%--------------------------------------------------------------------------
-estimated_params;
-BETTA;
-B;
-H;
-PHI2;
-RRA;
-KAPAtwo;
-ALFA;
-RHOR;
-BETTAPAI;
-BETTAY;
-MYYPS;
-MYZ;
-RHOA;
-RHOG;
-PAI;
-CONSxhr40;
-GoY;
-stderr eps_a;
-stderr eps_g;
-stderr eps_d;
-end;
-
-estimated_params_init(use_calibration);
-end;
-
-%--------------------------------------------------------------------------
-% Compare whether toolbox yields equivalent moments at second order
-%--------------------------------------------------------------------------
-% Note that we compare results for orderApp=1|2 and not for orderApp=3, because
-% there is a small error in the replication files of the original article in the 
-% computation of the covariance matrix of the extended innovations vector.
-% The authors have been contacted, fixed it, and report that the results 
-% change only slightly at orderApp=3 to what they report in the paper. At
-% orderApp=2 all is correct and so the following part tests whether we get 
-% the same model moments at the calibrated parameters (we do not optimize).
-% We compare it to the replication file RunGMM_Feedback_estim_RRA.m with the
-% following settings: orderApp=1|2, seOn=0, q_lag=10, weighting=1;
-% scaled=0; optimizer=0; estimator=1; momentSet=2;
-%
-% Output of the replication files for orderApp=1
-AndreasenEtAl.Q1 = 201778.9697;
-AndreasenEtAl.moments1 =[ % note that we reshuffeled to be compatible with our matched moments block
-    {[ 1]}    {'Ex'  }    {'Gr_C   '}    {'     '  }    {'0.024388'   }    {'0.023654'   }
-    {[ 2]}    {'Ex'  }    {'Gr_I   '}    {'     '  }    {'0.031046'   }    {'0.027719'   }
-    {[ 3]}    {'Ex'  }    {'Infl  ' }    {'     '  }    {'0.03757'    }    {'0.047415'   }
-    {[ 4]}    {'Ex'  }    {'r1    ' }    {'     '  }    {'0.056048'   }    {'0.083059'   }
-    {[ 5]}    {'Ex'  }    {'r40   ' }    {'     '  }    {'0.069929'   }    {'0.083059'   }
-    {[ 6]}    {'Ex'  }    {'xhr40  '}    {'     '  }    {'0.017237'   }    {'0'          }
-    {[ 7]}    {'Ex'  }    {'GoY    '}    {'     '  }    {'-1.5745'    }    {'-1.5745'    }
-    {[ 8]}    {'Ex'  }    {'hours  '}    {'     '  }    {'-0.043353'  }    {'-0.043245'  }
-    {[ 9]}    {'Exx' }    {'Gr_C   '}    {'Gr_C   '}    {'0.0013159'  }    {'0.0012253'  }
-    {[17]}    {'Exx' }    {'Gr_C   '}    {'Gr_I   '}    {'0.0021789'  }    {'0.0015117'  }
-    {[18]}    {'Exx' }    {'Gr_C   '}    {'Infl  ' }    {'0.00067495' }    {'0.00080078' }
-    {[19]}    {'Exx' }    {'Gr_C   '}    {'r1    ' }    {'0.0011655'  }    {'0.00182'    }
-    {[20]}    {'Exx' }    {'Gr_C   '}    {'r40   ' }    {'0.0015906'  }    {'0.001913'   }
-    {[21]}    {'Exx' }    {'Gr_C   '}    {'xhr40  '}    {'0.0020911'  }    {'0.0016326'  }
-    {[10]}    {'Exx' }    {'Gr_I   '}    {'Gr_I   '}    {'0.0089104'  }    {'0.0040112'  }
-    {[22]}    {'Exx' }    {'Gr_I   '}    {'Infl  ' }    {'0.00063139' }    {'0.00060604' }
-    {[23]}    {'Exx' }    {'Gr_I   '}    {'r1    ' }    {'0.0011031'  }    {'0.0021426'  }
-    {[24]}    {'Exx' }    {'Gr_I   '}    {'r40   ' }    {'0.0018445'  }    {'0.0022348'  }
-    {[25]}    {'Exx' }    {'Gr_I   '}    {'xhr40  '}    {'0.00095556' }    {'0.0039852'  }
-    {[11]}    {'Exx' }    {'Infl  ' }    {'Infl  ' }    {'0.0020268'  }    {'0.0030058'  }
-    {[26]}    {'Exx' }    {'Infl  ' }    {'r1    ' }    {'0.0025263'  }    {'0.0044951'  }
-    {[27]}    {'Exx' }    {'Infl  ' }    {'r40   ' }    {'0.0029126'  }    {'0.0042225'  }
-    {[28]}    {'Exx' }    {'Infl  ' }    {'xhr40  '}    {'-0.00077101'}    {'-0.0021222' }
-    {[12]}    {'Exx' }    {'r1    ' }    {'r1    ' }    {'0.0038708'  }    {'0.0074776'  }
-    {[29]}    {'Exx' }    {'r1    ' }    {'r40   ' }    {'0.0044773'  }    {'0.0071906'  }
-    {[30]}    {'Exx' }    {'r1    ' }    {'xhr40  '}    {'-0.00048202'}    {'-0.0006736' }
-    {[13]}    {'Exx' }    {'r40   ' }    {'r40   ' }    {'0.0054664'  }    {'0.0070599'  }
-    {[31]}    {'Exx' }    {'r40   ' }    {'xhr40  '}    {'0.00053864' }    {'-0.00036735'}
-    {[14]}    {'Exx' }    {'xhr40  '}    {'xhr40  '}    {'0.053097'   }    {'0.014516'   }
-    {[15]}    {'Exx' }    {'GoY    '}    {'GoY    '}    {'2.4863'     }    {'2.4866'     }
-    {[16]}    {'Exx' }    {'hours  '}    {'hours  '}    {'0.0018799'  }    {'0.0018713'  }
-    {[32]}    {'Exx1'}    {'Gr_C   '}    {'Gr_C   '}    {'0.00077917' }    {'0.00076856' }
-    {[33]}    {'Exx1'}    {'Gr_I   '}    {'Gr_I   '}    {'0.0050104'  }    {'0.002163'   }
-    {[34]}    {'Exx1'}    {'Infl  ' }    {'Infl  ' }    {'0.0019503'  }    {'0.0028078'  }
-    {[35]}    {'Exx1'}    {'r1    ' }    {'r1    ' }    {'0.0038509'  }    {'0.0074583'  }
-    {[36]}    {'Exx1'}    {'r40   ' }    {'r40   ' }    {'0.0054699'  }    {'0.0070551'  }
-    {[37]}    {'Exx1'}    {'xhr40  '}    {'xhr40  '}    {'-0.00098295'}    {'7.2164e-16' }
-    {[38]}    {'Exx1'}    {'GoY    '}    {'GoY    '}    {'2.4868'     }    {'2.4856'     }
-    {[39]}    {'Exx1'}    {'hours  '}    {'hours  '}    {'0.0018799'  }    {'0.0018708'  }
-];
-
-% Output of the replication files for orderApp=2
-AndreasenEtAl.Q2 = 59.3323;
-AndreasenEtAl.moments2 =[ % note that we reshuffeled to be compatible with our matched moments block
-    {[ 1]}    {'Ex'  }    {'Gr_C   '}    {'     '  }    {'0.024388'   }    {'0.023654'  }
-    {[ 2]}    {'Ex'  }    {'Gr_I   '}    {'     '  }    {'0.031046'   }    {'0.027719'  }
-    {[ 3]}    {'Ex'  }    {'Infl  ' }    {'     '  }    {'0.03757'    }    {'0.034565'  }
-    {[ 4]}    {'Ex'  }    {'r1    ' }    {'     '  }    {'0.056048'   }    {'0.056419'  }
-    {[ 5]}    {'Ex'  }    {'r40   ' }    {'     '  }    {'0.069929'   }    {'0.07087'   }
-    {[ 6]}    {'Ex'  }    {'xhr40  '}    {'     '  }    {'0.017237'   }    {'0.01517'   }
-    {[ 7]}    {'Ex'  }    {'GoY    '}    {'     '  }    {'-1.5745'    }    {'-1.5743'   }
-    {[ 8]}    {'Ex'  }    {'hours  '}    {'     '  }    {'-0.043353'  }    {'-0.043352' }
-    {[ 9]}    {'Exx' }    {'Gr_C   '}    {'Gr_C   '}    {'0.0013159'  }    {'0.0012464' }
-    {[17]}    {'Exx' }    {'Gr_C   '}    {'Gr_I   '}    {'0.0021789'  }    {'0.0015247' }
-    {[18]}    {'Exx' }    {'Gr_C   '}    {'Infl  ' }    {'0.00067495' }    {'0.0004867' }
-    {[19]}    {'Exx' }    {'Gr_C   '}    {'r1    ' }    {'0.0011655'  }    {'0.0011867' }
-    {[20]}    {'Exx' }    {'Gr_C   '}    {'r40   ' }    {'0.0015906'  }    {'0.0016146' }
-    {[21]}    {'Exx' }    {'Gr_C   '}    {'xhr40  '}    {'0.0020911'  }    {'0.0021395' }
-    {[10]}    {'Exx' }    {'Gr_I   '}    {'Gr_I   '}    {'0.0089104'  }    {'0.0043272' }
-    {[22]}    {'Exx' }    {'Gr_I   '}    {'Infl  ' }    {'0.00063139' }    {'0.00021752'}
-    {[23]}    {'Exx' }    {'Gr_I   '}    {'r1    ' }    {'0.0011031'  }    {'0.0013919' }
-    {[24]}    {'Exx' }    {'Gr_I   '}    {'r40   ' }    {'0.0018445'  }    {'0.0018899' }
-    {[25]}    {'Exx' }    {'Gr_I   '}    {'xhr40  '}    {'0.00095556' }    {'0.0037854' }
-    {[11]}    {'Exx' }    {'Infl  ' }    {'Infl  ' }    {'0.0020268'  }    {'0.0021043' }
-    {[26]}    {'Exx' }    {'Infl  ' }    {'r1    ' }    {'0.0025263'  }    {'0.0026571' }
-    {[27]}    {'Exx' }    {'Infl  ' }    {'r40   ' }    {'0.0029126'  }    {'0.0028566' }
-    {[28]}    {'Exx' }    {'Infl  ' }    {'xhr40  '}    {'-0.00077101'}    {'-0.0016279'}
-    {[12]}    {'Exx' }    {'r1    ' }    {'r1    ' }    {'0.0038708'  }    {'0.0039136' }
-    {[29]}    {'Exx' }    {'r1    ' }    {'r40   ' }    {'0.0044773'  }    {'0.0044118' }
-    {[30]}    {'Exx' }    {'r1    ' }    {'xhr40  '}    {'-0.00048202'}    {'0.00016791'}
-    {[13]}    {'Exx' }    {'r40   ' }    {'r40   ' }    {'0.0054664'  }    {'0.0052851' }
-    {[31]}    {'Exx' }    {'r40   ' }    {'xhr40  '}    {'0.00053864' }    {'0.00062143'}
-    {[14]}    {'Exx' }    {'xhr40  '}    {'xhr40  '}    {'0.053097'   }    {'0.018126'  }
-    {[15]}    {'Exx' }    {'GoY    '}    {'GoY    '}    {'2.4863'     }    {'2.4863'    }
-    {[16]}    {'Exx' }    {'hours  '}    {'hours  '}    {'0.0018799'  }    {'0.0018806' }
-    {[32]}    {'Exx1'}    {'Gr_C   '}    {'Gr_C   '}    {'0.00077917' }    {'0.00078586'}
-    {[33]}    {'Exx1'}    {'Gr_I   '}    {'Gr_I   '}    {'0.0050104'  }    {'0.0021519' }
-    {[34]}    {'Exx1'}    {'Infl  ' }    {'Infl  ' }    {'0.0019503'  }    {'0.0019046' }
-    {[35]}    {'Exx1'}    {'r1    ' }    {'r1    ' }    {'0.0038509'  }    {'0.0038939' }
-    {[36]}    {'Exx1'}    {'r40   ' }    {'r40   ' }    {'0.0054699'  }    {'0.0052792' }
-    {[37]}    {'Exx1'}    {'xhr40  '}    {'xhr40  '}    {'-0.00098295'}    {'0.00023012'}
-    {[38]}    {'Exx1'}    {'GoY    '}    {'GoY    '}    {'2.4868'     }    {'2.4852'    }
-    {[39]}    {'Exx1'}    {'hours  '}    {'hours  '}    {'0.0018799'  }    {'0.0018801' }
-];
-
-@#for orderApp in 1:2
-
-method_of_moments(
-          mom_method = GMM                    % method of moments method; possible values: GMM|SMM
-        , datafile   = 'AFVRR_data.mat'       % name of filename with data
-        , bartlett_kernel_lag = 10          % bandwith in optimal weighting matrix
-        , order = @{orderApp}                 % order of Taylor approximation in perturbation
-        , pruning                             % use pruned state space system at higher-order
-        % , verbose                           % display and store intermediate estimation results
-        , weighting_matrix = ['DIAGONAL']      % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename        
-        % , TeX                               % print TeX tables and graphics
-        % Optimization options that can be set by the user in the mod file, otherwise default values are provided        
-        %, huge_number=1D10                   % value for replacing the infinite bounds on parameters by finite numbers. Used by some optimizers for numerical reasons
-        , mode_compute = 0                    % specifies the optimizer for minimization of moments distance, note that by default there is a new optimizer        
-        , optim = ('TolFun', 1e-6
-                   ,'TolX', 1e-6
-                   ,'MaxIter', 3000
-                   ,'MaxFunEvals', 1D6
-                   ,'UseParallel' , 1
-                   %,'Jacobian' , 'on'
-                  )    % a list of NAME and VALUE pairs to set options for the optimization routines. Available options depend on mode_compute
-        %, silent_optimizer                  % run minimization of moments distance silently without displaying results or saving files in between
-        %, analytic_standard_errors
-        , se_tolx=1e-10
-);
-
-% Check results
-
-fprintf('****************************************************************\n')
-fprintf('Compare Results for perturbation order @{orderApp}\n')
-fprintf('****************************************************************\n')
-dev_Q            = AndreasenEtAl.Q@{orderApp} - oo_.mom.Q;
-dev_datamoments  = str2double(AndreasenEtAl.moments@{orderApp}(:,5)) - oo_.mom.data_moments;
-dev_modelmoments = str2double(AndreasenEtAl.moments@{orderApp}(:,6)) - oo_.mom.model_moments;
-
-% There is no table command in Octave
-% The table command also crashes on MATLAB R2014a because it does not like variable names
-if ~isoctave && ~matlab_ver_less_than('8.4')
-table([AndreasenEtAl.Q@{orderApp} ; str2double(AndreasenEtAl.moments@{orderApp}(:,5)) ; str2double(AndreasenEtAl.moments@{orderApp}(:,6))],...
-      [oo_.mom.Q                  ; oo_.mom.data_moments                              ; oo_.mom.model_moments                            ],...
-      [dev_Q                      ; dev_datamoments                                   ; dev_modelmoments                                 ],...
-      'VariableNames', {'Andreasen et al', 'Dynare', 'dev'},...
-      'RowNames', ['Q'; strcat('Data_', M_.matched_moments(:,4)); strcat('Model_', M_.matched_moments(:,4))])
-end
-
-if norm(dev_modelmoments)> 1e-4
-    warning('Something wrong in the computation of moments at order @{orderApp}')
-end
-
-@#endfor
-
-%--------------------------------------------------------------------------
-% Replicate estimation at orderApp=3
-%--------------------------------------------------------------------------
-@#ifdef DoEstimation
-method_of_moments(
-          mom_method = GMM                    % method of moments method; possible values: GMM|SMM
-        , datafile   = 'AFVRR_data.mat'       % name of filename with data
-        , bartlett_kernel_lag = 10            % bandwith in optimal weighting matrix
-        , order = 3                           % order of Taylor approximation in perturbation
-        , pruning                             % use pruned state space system at higher-order
-        % , verbose                           % display and store intermediate estimation results
-        , weighting_matrix = ['DIAGONAL', 'Optimal']      % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename        
-        % , TeX                               % print TeX tables and graphics
-        % Optimization options that can be set by the user in the mod file, otherwise default values are provided        
-        %, huge_number=1D10                   % value for replacing the infinite bounds on parameters by finite numbers. Used by some optimizers for numerical reasons
-        , mode_compute = 13                   % specifies the optimizer for minimization of moments distance, note that by default there is a new optimizer
-        , additional_optimizer_steps = [13]
-        , optim = ('TolFun', 1e-6
-                   ,'TolX', 1e-6
-                   ,'MaxIter', 3000
-                   ,'MaxFunEvals', 1D6
-                   ,'UseParallel' , 1
-                   %,'Jacobian' , 'on'
-                  )    % a list of NAME and VALUE pairs to set options for the optimization routines. Available options depend on mode_compute
-        %, silent_optimizer                  % run minimization of moments distance silently without displaying results or saving files in between
-        %, analytic_standard_errors
-        , se_tolx=1e-10
-);
-@#endif
+% DSGE model based on replication files of
+% Andreasen, Fernandez-Villaverde, Rubio-Ramirez (2018), The Pruned State-Space System for Non-Linear DSGE Models: Theory and Empirical Applications, Review of Economic Studies, 85, p. 1-49
+% Adapted for Dynare by Willi Mutschler (@wmutschl, willi@mutschler.eu), Jan 2021
+% =========================================================================
+% Copyright (C) 2021 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 <https://www.gnu.org/licenses/>.
+% =========================================================================
+
+% This is the model with Feedback M_FB
+% Original code RunGMM_Feedback_estim_RRA.m by Martin M. Andreasen, Jan 2016
+
+@#include "AFVRR_common.inc"
+
+%--------------------------------------------------------------------------
+% Parameter calibration taken from RunGMM_Feedback_estim_RRA.m
+%--------------------------------------------------------------------------
+% fixed parameters
+INHABIT   = 1;
+PHI1      = 4;
+PHI4      = 1;
+KAPAone   = 0;
+DELTA     = 0.025;
+THETA     = 0.36;
+ETA       = 6;
+CHI       = 0;
+BETTAxhr  = 0;
+BETTAxhr40= 0;
+RHOD      = 0;
+GAMA      = 0.9999;
+CONSxhr20 = 0;
+
+% estimated parameters
+BETTA     = 0.997007023687000;
+B         = 0.692501768577000;
+H         = 0.339214495653000;
+PHI2      = 0.688555040951000;
+RRA       = 24.346514272871001;
+KAPAtwo   = 10.018421876923000;
+ALFA      = 0.792507553312000;
+RHOR      = 0.849194030384000;
+BETTAPAI  = 2.060579322980000;
+BETTAY    = 0.220573712342000;
+MYYPS     = 1.001016690133000;
+MYZ       = 1.005356313981000;
+RHOA      = 0.784141391843000;
+RHOG      = 0.816924540497000;
+PAI       = 1.011924196487000;
+CONSxhr40 = 0.878774662208000;
+GoY       = 0.207110300602000;
+STDA      = 0.013024450606000;
+STDG      = 0.051049871928000;
+STDD      = 0.008877423780000;
+
+% endogenous parameters set via steady state, no need to initialize
+%PHIzero  = ;
+%AA       = ;
+%PHI3     = ;
+%negVf    = ;
+
+model_diagnostics;
+% Model diagnostics show that some parameters are endogenously determined
+% via the steady state, so we run steady to calibrate all parameters
+steady;
+model_diagnostics;
+% Now all parameters are determined
+
+resid;
+check;
+
+%--------------------------------------------------------------------------
+% Shock distribution
+%--------------------------------------------------------------------------
+shocks;
+var eps_a = STDA^2;
+var eps_d = STDD^2;
+var eps_g = STDG^2;
+end;
+
+%--------------------------------------------------------------------------
+% Estimated Params block - these parameters will be estimated, we
+% initialize at calibrated values
+%--------------------------------------------------------------------------
+estimated_params;
+BETTA;
+B;
+H;
+PHI2;
+RRA;
+KAPAtwo;
+ALFA;
+RHOR;
+BETTAPAI;
+BETTAY;
+MYYPS;
+MYZ;
+RHOA;
+RHOG;
+PAI;
+CONSxhr40;
+GoY;
+stderr eps_a;
+stderr eps_g;
+stderr eps_d;
+end;
+
+estimated_params_init(use_calibration);
+end;
+
+%--------------------------------------------------------------------------
+% Compare whether toolbox yields equivalent moments at second order
+%--------------------------------------------------------------------------
+% Note that we compare results for orderApp=1|2 and not for orderApp=3, because
+% there is a small error in the replication files of the original article in the 
+% computation of the covariance matrix of the extended innovations vector.
+% The authors have been contacted, fixed it, and report that the results 
+% change only slightly at orderApp=3 to what they report in the paper. At
+% orderApp=2 all is correct and so the following part tests whether we get 
+% the same model moments at the calibrated parameters (we do not optimize).
+% We compare it to the replication file RunGMM_Feedback_estim_RRA.m with the
+% following settings: orderApp=1|2, seOn=0, q_lag=10, weighting=1;
+% scaled=0; optimizer=0; estimator=1; momentSet=2;
+%
+% Output of the replication files for orderApp=1
+AndreasenEtAl.Q1 = 201778.9697;
+AndreasenEtAl.moments1 =[ % note that we reshuffeled to be compatible with our matched moments block
+    {[ 1]}    {'Ex'  }    {'Gr_C   '}    {'     '  }    {'0.024388'   }    {'0.023654'   }
+    {[ 2]}    {'Ex'  }    {'Gr_I   '}    {'     '  }    {'0.031046'   }    {'0.027719'   }
+    {[ 3]}    {'Ex'  }    {'Infl  ' }    {'     '  }    {'0.03757'    }    {'0.047415'   }
+    {[ 4]}    {'Ex'  }    {'r1    ' }    {'     '  }    {'0.056048'   }    {'0.083059'   }
+    {[ 5]}    {'Ex'  }    {'r40   ' }    {'     '  }    {'0.069929'   }    {'0.083059'   }
+    {[ 6]}    {'Ex'  }    {'xhr40  '}    {'     '  }    {'0.017237'   }    {'0'          }
+    {[ 7]}    {'Ex'  }    {'GoY    '}    {'     '  }    {'-1.5745'    }    {'-1.5745'    }
+    {[ 8]}    {'Ex'  }    {'hours  '}    {'     '  }    {'-0.043353'  }    {'-0.043245'  }
+    {[ 9]}    {'Exx' }    {'Gr_C   '}    {'Gr_C   '}    {'0.0013159'  }    {'0.0012253'  }
+    {[10]}    {'Exx' }    {'Gr_C   '}    {'Gr_I   '}    {'0.0021789'  }    {'0.0015117'  }
+    {[11]}    {'Exx' }    {'Gr_C   '}    {'Infl  ' }    {'0.00067495' }    {'0.00080078' }
+    {[12]}    {'Exx' }    {'Gr_C   '}    {'r1    ' }    {'0.0011655'  }    {'0.00182'    }
+    {[13]}    {'Exx' }    {'Gr_C   '}    {'r40   ' }    {'0.0015906'  }    {'0.001913'   }
+    {[14]}    {'Exx' }    {'Gr_C   '}    {'xhr40  '}    {'0.0020911'  }    {'0.0016326'  }
+    {[15]}    {'Exx' }    {'Gr_I   '}    {'Gr_I   '}    {'0.0089104'  }    {'0.0040112'  }
+    {[16]}    {'Exx' }    {'Gr_I   '}    {'Infl  ' }    {'0.00063139' }    {'0.00060604' }
+    {[17]}    {'Exx' }    {'Gr_I   '}    {'r1    ' }    {'0.0011031'  }    {'0.0021426'  }
+    {[18]}    {'Exx' }    {'Gr_I   '}    {'r40   ' }    {'0.0018445'  }    {'0.0022348'  }
+    {[19]}    {'Exx' }    {'Gr_I   '}    {'xhr40  '}    {'0.00095556' }    {'0.0039852'  }
+    {[20]}    {'Exx' }    {'Infl  ' }    {'Infl  ' }    {'0.0020268'  }    {'0.0030058'  }
+    {[21]}    {'Exx' }    {'Infl  ' }    {'r1    ' }    {'0.0025263'  }    {'0.0044951'  }
+    {[22]}    {'Exx' }    {'Infl  ' }    {'r40   ' }    {'0.0029126'  }    {'0.0042225'  }
+    {[23]}    {'Exx' }    {'Infl  ' }    {'xhr40  '}    {'-0.00077101'}    {'-0.0021222' }
+    {[24]}    {'Exx' }    {'r1    ' }    {'r1    ' }    {'0.0038708'  }    {'0.0074776'  }
+    {[25]}    {'Exx' }    {'r1    ' }    {'r40   ' }    {'0.0044773'  }    {'0.0071906'  }
+    {[26]}    {'Exx' }    {'r1    ' }    {'xhr40  '}    {'-0.00048202'}    {'-0.0006736' }
+    {[27]}    {'Exx' }    {'r40   ' }    {'r40   ' }    {'0.0054664'  }    {'0.0070599'  }
+    {[28]}    {'Exx' }    {'r40   ' }    {'xhr40  '}    {'0.00053864' }    {'-0.00036735'}
+    {[29]}    {'Exx' }    {'xhr40  '}    {'xhr40  '}    {'0.053097'   }    {'0.014516'   }
+    {[30]}    {'Exx' }    {'GoY    '}    {'GoY    '}    {'2.4863'     }    {'2.4866'     }
+    {[31]}    {'Exx' }    {'hours  '}    {'hours  '}    {'0.0018799'  }    {'0.0018713'  }
+    {[32]}    {'Exx1'}    {'Gr_C   '}    {'Gr_C   '}    {'0.00077917' }    {'0.00076856' }
+    {[33]}    {'Exx1'}    {'Gr_I   '}    {'Gr_I   '}    {'0.0050104'  }    {'0.002163'   }
+    {[34]}    {'Exx1'}    {'Infl  ' }    {'Infl  ' }    {'0.0019503'  }    {'0.0028078'  }
+    {[35]}    {'Exx1'}    {'r1    ' }    {'r1    ' }    {'0.0038509'  }    {'0.0074583'  }
+    {[36]}    {'Exx1'}    {'r40   ' }    {'r40   ' }    {'0.0054699'  }    {'0.0070551'  }
+    {[37]}    {'Exx1'}    {'xhr40  '}    {'xhr40  '}    {'-0.00098295'}    {'7.2164e-16' }
+    {[38]}    {'Exx1'}    {'GoY    '}    {'GoY    '}    {'2.4868'     }    {'2.4856'     }
+    {[39]}    {'Exx1'}    {'hours  '}    {'hours  '}    {'0.0018799'  }    {'0.0018708'  }
+];
+
+% Output of the replication files for orderApp=2
+AndreasenEtAl.Q2 = 59.3323;
+AndreasenEtAl.moments2 =[ % note that we reshuffeled to be compatible with our matched moments block
+    {[ 1]}    {'Ex'  }    {'Gr_C   '}    {'     '  }    {'0.024388'   }    {'0.023654'  }
+    {[ 2]}    {'Ex'  }    {'Gr_I   '}    {'     '  }    {'0.031046'   }    {'0.027719'  }
+    {[ 3]}    {'Ex'  }    {'Infl  ' }    {'     '  }    {'0.03757'    }    {'0.034565'  }
+    {[ 4]}    {'Ex'  }    {'r1    ' }    {'     '  }    {'0.056048'   }    {'0.056419'  }
+    {[ 5]}    {'Ex'  }    {'r40   ' }    {'     '  }    {'0.069929'   }    {'0.07087'   }
+    {[ 6]}    {'Ex'  }    {'xhr40  '}    {'     '  }    {'0.017237'   }    {'0.01517'   }
+    {[ 7]}    {'Ex'  }    {'GoY    '}    {'     '  }    {'-1.5745'    }    {'-1.5743'   }
+    {[ 8]}    {'Ex'  }    {'hours  '}    {'     '  }    {'-0.043353'  }    {'-0.043352' }
+    {[ 9]}    {'Exx' }    {'Gr_C   '}    {'Gr_C   '}    {'0.0013159'  }    {'0.0012464' }
+    {[10]}    {'Exx' }    {'Gr_C   '}    {'Gr_I   '}    {'0.0021789'  }    {'0.0015247' }
+    {[11]}    {'Exx' }    {'Gr_C   '}    {'Infl  ' }    {'0.00067495' }    {'0.0004867' }
+    {[12]}    {'Exx' }    {'Gr_C   '}    {'r1    ' }    {'0.0011655'  }    {'0.0011867' }
+    {[13]}    {'Exx' }    {'Gr_C   '}    {'r40   ' }    {'0.0015906'  }    {'0.0016146' }
+    {[14]}    {'Exx' }    {'Gr_C   '}    {'xhr40  '}    {'0.0020911'  }    {'0.0021395' }
+    {[15]}    {'Exx' }    {'Gr_I   '}    {'Gr_I   '}    {'0.0089104'  }    {'0.0043272' }
+    {[16]}    {'Exx' }    {'Gr_I   '}    {'Infl  ' }    {'0.00063139' }    {'0.00021752'}
+    {[17]}    {'Exx' }    {'Gr_I   '}    {'r1    ' }    {'0.0011031'  }    {'0.0013919' }
+    {[18]}    {'Exx' }    {'Gr_I   '}    {'r40   ' }    {'0.0018445'  }    {'0.0018899' }
+    {[19]}    {'Exx' }    {'Gr_I   '}    {'xhr40  '}    {'0.00095556' }    {'0.0037854' }
+    {[20]}    {'Exx' }    {'Infl  ' }    {'Infl  ' }    {'0.0020268'  }    {'0.0021043' }
+    {[21]}    {'Exx' }    {'Infl  ' }    {'r1    ' }    {'0.0025263'  }    {'0.0026571' }
+    {[22]}    {'Exx' }    {'Infl  ' }    {'r40   ' }    {'0.0029126'  }    {'0.0028566' }
+    {[23]}    {'Exx' }    {'Infl  ' }    {'xhr40  '}    {'-0.00077101'}    {'-0.0016279'}
+    {[24]}    {'Exx' }    {'r1    ' }    {'r1    ' }    {'0.0038708'  }    {'0.0039136' }
+    {[25]}    {'Exx' }    {'r1    ' }    {'r40   ' }    {'0.0044773'  }    {'0.0044118' }
+    {[26]}    {'Exx' }    {'r1    ' }    {'xhr40  '}    {'-0.00048202'}    {'0.00016791'}
+    {[27]}    {'Exx' }    {'r40   ' }    {'r40   ' }    {'0.0054664'  }    {'0.0052851' }
+    {[28]}    {'Exx' }    {'r40   ' }    {'xhr40  '}    {'0.00053864' }    {'0.00062143'}
+    {[29]}    {'Exx' }    {'xhr40  '}    {'xhr40  '}    {'0.053097'   }    {'0.018126'  }
+    {[30]}    {'Exx' }    {'GoY    '}    {'GoY    '}    {'2.4863'     }    {'2.4863'    }
+    {[31]}    {'Exx' }    {'hours  '}    {'hours  '}    {'0.0018799'  }    {'0.0018806' }
+    {[32]}    {'Exx1'}    {'Gr_C   '}    {'Gr_C   '}    {'0.00077917' }    {'0.00078586'}
+    {[33]}    {'Exx1'}    {'Gr_I   '}    {'Gr_I   '}    {'0.0050104'  }    {'0.0021519' }
+    {[34]}    {'Exx1'}    {'Infl  ' }    {'Infl  ' }    {'0.0019503'  }    {'0.0019046' }
+    {[35]}    {'Exx1'}    {'r1    ' }    {'r1    ' }    {'0.0038509'  }    {'0.0038939' }
+    {[36]}    {'Exx1'}    {'r40   ' }    {'r40   ' }    {'0.0054699'  }    {'0.0052792' }
+    {[37]}    {'Exx1'}    {'xhr40  '}    {'xhr40  '}    {'-0.00098295'}    {'0.00023012'}
+    {[38]}    {'Exx1'}    {'GoY    '}    {'GoY    '}    {'2.4868'     }    {'2.4852'    }
+    {[39]}    {'Exx1'}    {'hours  '}    {'hours  '}    {'0.0018799'  }    {'0.0018801' }
+];
+
+@#for orderApp in 1:2
+
+method_of_moments(
+          mom_method = GMM                    % method of moments method; possible values: GMM|SMM
+        , datafile   = 'AFVRR_data.mat'       % name of filename with data
+        , bartlett_kernel_lag = 10          % bandwith in optimal weighting matrix
+        , order = @{orderApp}                 % order of Taylor approximation in perturbation
+        , pruning                             % use pruned state space system at higher-order
+        % , verbose                           % display and store intermediate estimation results
+        , weighting_matrix = ['DIAGONAL']      % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename        
+        % , TeX                               % print TeX tables and graphics
+        % Optimization options that can be set by the user in the mod file, otherwise default values are provided        
+        %, huge_number=1D10                   % value for replacing the infinite bounds on parameters by finite numbers. Used by some optimizers for numerical reasons
+        , mode_compute = 0                    % specifies the optimizer for minimization of moments distance, note that by default there is a new optimizer        
+        , optim = ('TolFun', 1e-6
+                   ,'TolX', 1e-6
+                   ,'MaxIter', 3000
+                   ,'MaxFunEvals', 1D6
+                   ,'UseParallel' , 1
+                   %,'Jacobian' , 'on'
+                  )    % a list of NAME and VALUE pairs to set options for the optimization routines. Available options depend on mode_compute
+        %, silent_optimizer                  % run minimization of moments distance silently without displaying results or saving files in between
+        %, analytic_standard_errors
+        , se_tolx=1e-10
+);
+
+% Check results
+
+fprintf('****************************************************************\n')
+fprintf('Compare Results for perturbation order @{orderApp}\n')
+fprintf('****************************************************************\n')
+dev_Q            = AndreasenEtAl.Q@{orderApp} - oo_.mom.Q;
+dev_datamoments  = str2double(AndreasenEtAl.moments@{orderApp}(:,5)) - oo_.mom.data_moments;
+dev_modelmoments = str2double(AndreasenEtAl.moments@{orderApp}(:,6)) - oo_.mom.model_moments;
+
+% There is no table command in Octave
+% The table command also crashes on MATLAB R2014a because it does not like variable names
+if ~isoctave && ~matlab_ver_less_than('8.4')
+table([AndreasenEtAl.Q@{orderApp} ; str2double(AndreasenEtAl.moments@{orderApp}(:,5)) ; str2double(AndreasenEtAl.moments@{orderApp}(:,6))],...
+      [oo_.mom.Q                  ; oo_.mom.data_moments                              ; oo_.mom.model_moments                            ],...
+      [dev_Q                      ; dev_datamoments                                   ; dev_modelmoments                                 ],...
+      'VariableNames', {'Andreasen et al', 'Dynare', 'dev'})
+end
+
+if norm(dev_modelmoments)> 1e-4
+    warning('Something wrong in the computation of moments at order @{orderApp}')
+end
+
+@#endfor
+
+%--------------------------------------------------------------------------
+% Replicate estimation at orderApp=3
+%--------------------------------------------------------------------------
+@#ifdef DoEstimation
+method_of_moments(
+          mom_method = GMM                    % method of moments method; possible values: GMM|SMM
+        , datafile   = 'AFVRR_data.mat'       % name of filename with data
+        , bartlett_kernel_lag = 10            % bandwith in optimal weighting matrix
+        , order = 3                           % order of Taylor approximation in perturbation
+        , pruning                             % use pruned state space system at higher-order
+        % , verbose                           % display and store intermediate estimation results
+        , weighting_matrix = ['DIAGONAL', 'Optimal']      % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename        
+        % , TeX                               % print TeX tables and graphics
+        % Optimization options that can be set by the user in the mod file, otherwise default values are provided        
+        %, huge_number=1D10                   % value for replacing the infinite bounds on parameters by finite numbers. Used by some optimizers for numerical reasons
+        , mode_compute = 13                   % specifies the optimizer for minimization of moments distance, note that by default there is a new optimizer
+        , additional_optimizer_steps = [13]
+        , optim = ('TolFun', 1e-6
+                   ,'TolX', 1e-6
+                   ,'MaxIter', 3000
+                   ,'MaxFunEvals', 1D6
+                   ,'UseParallel' , 1
+                   %,'Jacobian' , 'on'
+                  )    % a list of NAME and VALUE pairs to set options for the optimization routines. Available options depend on mode_compute
+        %, silent_optimizer                  % run minimization of moments distance silently without displaying results or saving files in between
+        %, analytic_standard_errors
+        , se_tolx=1e-10
+);
+@#endif

tests/estimation/method_of_moments/AFVRR/AFVRR_MFB_RRA.mod

@@ -145,28 +145,28 @@ AndreasenEtAl.moments1 =[ % note that we reshuffeled to be compatible with our m
     {[ 7]}    {'Ex'  }    {'GoY    '}    {'     '  }    {'-1.5745'    }    {'-1.5746'    }
     {[ 8]}    {'Ex'  }    {'hours  '}    {'     '  }    {'-0.043353'  }    {'-0.043299'  }
     {[ 9]}    {'Exx' }    {'Gr_C   '}    {'Gr_C   '}    {'0.0013159'  }    {'0.0012763'  }
-    {[17]}    {'Exx' }    {'Gr_C   '}    {'Gr_I   '}    {'0.0021789'  }    {'0.0017759'  }
-    {[18]}    {'Exx' }    {'Gr_C   '}    {'Infl  ' }    {'0.00067495' }    {'0.00077354' }
-    {[19]}    {'Exx' }    {'Gr_C   '}    {'r1    ' }    {'0.0011655'  }    {'0.0016538'  }
-    {[20]}    {'Exx' }    {'Gr_C   '}    {'r40   ' }    {'0.0015906'  }    {'0.0017949'  }
-    {[21]}    {'Exx' }    {'Gr_C   '}    {'xhr40  '}    {'0.0020911'  }    {'0.0017847'  }
-    {[10]}    {'Exx' }    {'Gr_I   '}    {'Gr_I   '}    {'0.0089104'  }    {'0.0053424'  }
-    {[22]}    {'Exx' }    {'Gr_I   '}    {'Infl  ' }    {'0.00063139' }    {'0.00064897' }
-    {[23]}    {'Exx' }    {'Gr_I   '}    {'r1    ' }    {'0.0011031'  }    {'0.0019533'  }
-    {[24]}    {'Exx' }    {'Gr_I   '}    {'r40   ' }    {'0.0018445'  }    {'0.0020602'  }
-    {[25]}    {'Exx' }    {'Gr_I   '}    {'xhr40  '}    {'0.00095556' }    {'0.0064856'  }
-    {[11]}    {'Exx' }    {'Infl  ' }    {'Infl  ' }    {'0.0020268'  }    {'0.0020922'  }
-    {[26]}    {'Exx' }    {'Infl  ' }    {'r1    ' }    {'0.0025263'  }    {'0.0036375'  }
-    {[27]}    {'Exx' }    {'Infl  ' }    {'r40   ' }    {'0.0029126'  }    {'0.0034139'  }
-    {[28]}    {'Exx' }    {'Infl  ' }    {'xhr40  '}    {'-0.00077101'}    {'-0.0011665' }
-    {[12]}    {'Exx' }    {'r1    ' }    {'r1    ' }    {'0.0038708'  }    {'0.0066074'  }
-    {[29]}    {'Exx' }    {'r1    ' }    {'r40   ' }    {'0.0044773'  }    {'0.0062959'  }
-    {[30]}    {'Exx' }    {'r1    ' }    {'xhr40  '}    {'-0.00048202'}    {'-0.00075499'}
-    {[13]}    {'Exx' }    {'r40   ' }    {'r40   ' }    {'0.0054664'  }    {'0.0061801'  }
-    {[31]}    {'Exx' }    {'r40   ' }    {'xhr40  '}    {'0.00053864' }    {'-0.00030456'}
-    {[14]}    {'Exx' }    {'xhr40  '}    {'xhr40  '}    {'0.053097'   }    {'0.012048'   }
-    {[15]}    {'Exx' }    {'GoY    '}    {'GoY    '}    {'2.4863'     }    {'2.4872'     }
-    {[16]}    {'Exx' }    {'hours  '}    {'hours  '}    {'0.0018799'  }    {'0.0018759'  }
+    {[10]}    {'Exx' }    {'Gr_C   '}    {'Gr_I   '}    {'0.0021789'  }    {'0.0017759'  }
+    {[11]}    {'Exx' }    {'Gr_C   '}    {'Infl  ' }    {'0.00067495' }    {'0.00077354' }
+    {[12]}    {'Exx' }    {'Gr_C   '}    {'r1    ' }    {'0.0011655'  }    {'0.0016538'  }
+    {[13]}    {'Exx' }    {'Gr_C   '}    {'r40   ' }    {'0.0015906'  }    {'0.0017949'  }
+    {[14]}    {'Exx' }    {'Gr_C   '}    {'xhr40  '}    {'0.0020911'  }    {'0.0017847'  }
+    {[15]}    {'Exx' }    {'Gr_I   '}    {'Gr_I   '}    {'0.0089104'  }    {'0.0053424'  }
+    {[16]}    {'Exx' }    {'Gr_I   '}    {'Infl  ' }    {'0.00063139' }    {'0.00064897' }
+    {[17]}    {'Exx' }    {'Gr_I   '}    {'r1    ' }    {'0.0011031'  }    {'0.0019533'  }
+    {[18]}    {'Exx' }    {'Gr_I   '}    {'r40   ' }    {'0.0018445'  }    {'0.0020602'  }
+    {[19]}    {'Exx' }    {'Gr_I   '}    {'xhr40  '}    {'0.00095556' }    {'0.0064856'  }
+    {[20]}    {'Exx' }    {'Infl  ' }    {'Infl  ' }    {'0.0020268'  }    {'0.0020922'  }
+    {[21]}    {'Exx' }    {'Infl  ' }    {'r1    ' }    {'0.0025263'  }    {'0.0036375'  }
+    {[22]}    {'Exx' }    {'Infl  ' }    {'r40   ' }    {'0.0029126'  }    {'0.0034139'  }
+    {[23]}    {'Exx' }    {'Infl  ' }    {'xhr40  '}    {'-0.00077101'}    {'-0.0011665' }
+    {[24]}    {'Exx' }    {'r1    ' }    {'r1    ' }    {'0.0038708'  }    {'0.0066074'  }
+    {[25]}    {'Exx' }    {'r1    ' }    {'r40   ' }    {'0.0044773'  }    {'0.0062959'  }
+    {[26]}    {'Exx' }    {'r1    ' }    {'xhr40  '}    {'-0.00048202'}    {'-0.00075499'}
+    {[27]}    {'Exx' }    {'r40   ' }    {'r40   ' }    {'0.0054664'  }    {'0.0061801'  }
+    {[28]}    {'Exx' }    {'r40   ' }    {'xhr40  '}    {'0.00053864' }    {'-0.00030456'}
+    {[29]}    {'Exx' }    {'xhr40  '}    {'xhr40  '}    {'0.053097'   }    {'0.012048'   }
+    {[30]}    {'Exx' }    {'GoY    '}    {'GoY    '}    {'2.4863'     }    {'2.4872'     }
+    {[31]}    {'Exx' }    {'hours  '}    {'hours  '}    {'0.0018799'  }    {'0.0018759'  }
     {[32]}    {'Exx1'}    {'Gr_C   '}    {'Gr_C   '}    {'0.00077917' }    {'0.00080528' }
     {[33]}    {'Exx1'}    {'Gr_I   '}    {'Gr_I   '}    {'0.0050104'  }    {'0.0017036'  }
     {[34]}    {'Exx1'}    {'Infl  ' }    {'Infl  ' }    {'0.0019503'  }    {'0.0020185'  }
@@ -189,28 +189,28 @@ AndreasenEtAl.moments2 =[ % note that we reshuffeled to be compatible with our m
     {[ 7]}    {'Ex'  }    {'GoY    '}    {'     '  }    {'-1.5745'    }    {'-1.574'     }
     {[ 8]}    {'Ex'  }    {'hours  '}    {'     '  }    {'-0.043353'  }    {'-0.043351'  }
     {[ 9]}    {'Exx' }    {'Gr_C   '}    {'Gr_C   '}    {'0.0013159'  }    {'0.0012917'  }
-    {[17]}    {'Exx' }    {'Gr_C   '}    {'Gr_I   '}    {'0.0021789'  }    {'0.0017862'  }
-    {[18]}    {'Exx' }    {'Gr_C   '}    {'Infl  ' }    {'0.00067495' }    {'0.00061078' }
-    {[19]}    {'Exx' }    {'Gr_C   '}    {'r1    ' }    {'0.0011655'  }    {'0.0011494'  }
-    {[20]}    {'Exx' }    {'Gr_C   '}    {'r40   ' }    {'0.0015906'  }    {'0.0016149'  }
-    {[21]}    {'Exx' }    {'Gr_C   '}    {'xhr40  '}    {'0.0020911'  }    {'0.002203'   }
-    {[10]}    {'Exx' }    {'Gr_I   '}    {'Gr_I   '}    {'0.0089104'  }    {'0.0054317'  }
-    {[22]}    {'Exx' }    {'Gr_I   '}    {'Infl  ' }    {'0.00063139' }    {'0.00045278' }
-    {[23]}    {'Exx' }    {'Gr_I   '}    {'r1    ' }    {'0.0011031'  }    {'0.0013672'  }
-    {[24]}    {'Exx' }    {'Gr_I   '}    {'r40   ' }    {'0.0018445'  }    {'0.0018557'  }
-    {[25]}    {'Exx' }    {'Gr_I   '}    {'xhr40  '}    {'0.00095556' }    {'0.0067742'  }
-    {[11]}    {'Exx' }    {'Infl  ' }    {'Infl  ' }    {'0.0020268'  }    {'0.0016583'  }
-    {[26]}    {'Exx' }    {'Infl  ' }    {'r1    ' }    {'0.0025263'  }    {'0.0024521'  }
-    {[27]}    {'Exx' }    {'Infl  ' }    {'r40   ' }    {'0.0029126'  }    {'0.002705'   }
-    {[28]}    {'Exx' }    {'Infl  ' }    {'xhr40  '}    {'-0.00077101'}    {'-0.00065007'}
-    {[12]}    {'Exx' }    {'r1    ' }    {'r1    ' }    {'0.0038708'  }    {'0.0038274'  }
-    {[29]}    {'Exx' }    {'r1    ' }    {'r40   ' }    {'0.0044773'  }    {'0.004297'   }
-    {[30]}    {'Exx' }    {'r1    ' }    {'xhr40  '}    {'-0.00048202'}    {'6.3243e-05' }
-    {[13]}    {'Exx' }    {'r40   ' }    {'r40   ' }    {'0.0054664'  }    {'0.0051686'  }
-    {[31]}    {'Exx' }    {'r40   ' }    {'xhr40  '}    {'0.00053864' }    {'0.00066645' }
-    {[14]}    {'Exx' }    {'xhr40  '}    {'xhr40  '}    {'0.053097'   }    {'0.013543'   }
-    {[15]}    {'Exx' }    {'GoY    '}    {'GoY    '}    {'2.4863'     }    {'2.4858'     }
-    {[16]}    {'Exx' }    {'hours  '}    {'hours  '}    {'0.0018799'  }    {'0.0018804'  }
+    {[10]}    {'Exx' }    {'Gr_C   '}    {'Gr_I   '}    {'0.0021789'  }    {'0.0017862'  }
+    {[11]}    {'Exx' }    {'Gr_C   '}    {'Infl  ' }    {'0.00067495' }    {'0.00061078' }
+    {[12]}    {'Exx' }    {'Gr_C   '}    {'r1    ' }    {'0.0011655'  }    {'0.0011494'  }
+    {[13]}    {'Exx' }    {'Gr_C   '}    {'r40   ' }    {'0.0015906'  }    {'0.0016149'  }
+    {[14]}    {'Exx' }    {'Gr_C   '}    {'xhr40  '}    {'0.0020911'  }    {'0.002203'   }
+    {[15]}    {'Exx' }    {'Gr_I   '}    {'Gr_I   '}    {'0.0089104'  }    {'0.0054317'  }
+    {[16]}    {'Exx' }    {'Gr_I   '}    {'Infl  ' }    {'0.00063139' }    {'0.00045278' }
+    {[17]}    {'Exx' }    {'Gr_I   '}    {'r1    ' }    {'0.0011031'  }    {'0.0013672'  }
+    {[18]}    {'Exx' }    {'Gr_I   '}    {'r40   ' }    {'0.0018445'  }    {'0.0018557'  }
+    {[19]}    {'Exx' }    {'Gr_I   '}    {'xhr40  '}    {'0.00095556' }    {'0.0067742'  }
+    {[20]}    {'Exx' }    {'Infl  ' }    {'Infl  ' }    {'0.0020268'  }    {'0.0016583'  }
+    {[21]}    {'Exx' }    {'Infl  ' }    {'r1    ' }    {'0.0025263'  }    {'0.0024521'  }
+    {[22]}    {'Exx' }    {'Infl  ' }    {'r40   ' }    {'0.0029126'  }    {'0.002705'   }
+    {[23]}    {'Exx' }    {'Infl  ' }    {'xhr40  '}    {'-0.00077101'}    {'-0.00065007'}
+    {[24]}    {'Exx' }    {'r1    ' }    {'r1    ' }    {'0.0038708'  }    {'0.0038274'  }
+    {[25]}    {'Exx' }    {'r1    ' }    {'r40   ' }    {'0.0044773'  }    {'0.004297'   }
+    {[26]}    {'Exx' }    {'r1    ' }    {'xhr40  '}    {'-0.00048202'}    {'6.3243e-05' }
+    {[27]}    {'Exx' }    {'r40   ' }    {'r40   ' }    {'0.0054664'  }    {'0.0051686'  }
+    {[28]}    {'Exx' }    {'r40   ' }    {'xhr40  '}    {'0.00053864' }    {'0.00066645' }
+    {[29]}    {'Exx' }    {'xhr40  '}    {'xhr40  '}    {'0.053097'   }    {'0.013543'   }
+    {[30]}    {'Exx' }    {'GoY    '}    {'GoY    '}    {'2.4863'     }    {'2.4858'     }
+    {[31]}    {'Exx' }    {'hours  '}    {'hours  '}    {'0.0018799'  }    {'0.0018804'  }
     {[32]}    {'Exx1'}    {'Gr_C   '}    {'Gr_C   '}    {'0.00077917' }    {'0.00081772' }
     {[33]}    {'Exx1'}    {'Gr_I   '}    {'Gr_I   '}    {'0.0050104'  }    {'0.0017106'  }
     {[34]}    {'Exx1'}    {'Infl  ' }    {'Infl  ' }    {'0.0019503'  }    {'0.0015835'  }
@@ -262,8 +262,7 @@ if ~isoctave && ~matlab_ver_less_than('8.4')
 table([AndreasenEtAl.Q@{orderApp} ; str2double(AndreasenEtAl.moments@{orderApp}(:,5)) ; str2double(AndreasenEtAl.moments@{orderApp}(:,6))],...
       [oo_.mom.Q                  ; oo_.mom.data_moments                              ; oo_.mom.model_moments                            ],...
       [dev_Q                      ; dev_datamoments                                   ; dev_modelmoments                                 ],...
-      'VariableNames', {'Andreasen et al', 'Dynare', 'dev'},...
-      'RowNames', ['Q'; strcat('Data_', M_.matched_moments(:,4)); strcat('Model_', M_.matched_moments(:,4))])
+      'VariableNames', {'Andreasen et al', 'Dynare', 'dev'})
 end
 
 if norm(dev_modelmoments)> 1e-4

tests/estimation/method_of_moments/AFVRR/AFVRR_common.inc

@@ -446,7 +446,7 @@ end;
 % These are the moments used in the paper; corresponds to momentSet=2 in the replication files
 
 matched_moments;
-%mean
+%first moments: all
 Obs_Gr_C;
 Obs_Gr_I;
 Obs_Infl;
@@ -456,53 +456,39 @@ Obs_xhr40;
 Obs_GoY;
 Obs_hours;
 
-% all variances
+%second moments
+% (i) all variances, (2) all covariances excluding GoY and hours, (3) own first autocovariances
 Obs_Gr_C*Obs_Gr_C;
-Obs_Gr_I*Obs_Gr_I;
-Obs_Infl*Obs_Infl;
-Obs_r1*Obs_r1;
-Obs_r40*Obs_r40;
-Obs_xhr40*Obs_xhr40;
-Obs_GoY*Obs_GoY;
-Obs_hours*Obs_hours;
-
-% covariance excluding GoY and hours
 Obs_Gr_C*Obs_Gr_I;
 Obs_Gr_C*Obs_Infl;
 Obs_Gr_C*Obs_r1;
 Obs_Gr_C*Obs_r40;
 Obs_Gr_C*Obs_xhr40;
-%Obs_Gr_C*Obs_GoY;
-%Obs_Gr_C*Obs_hours;
 
+Obs_Gr_I*Obs_Gr_I;
 Obs_Gr_I*Obs_Infl;
 Obs_Gr_I*Obs_r1;
 Obs_Gr_I*Obs_r40;
 Obs_Gr_I*Obs_xhr40;
-%Obs_Gr_I*Obs_GoY;
-%Obs_Gr_I*Obs_hours;
 
+Obs_Infl*Obs_Infl;
 Obs_Infl*Obs_r1;
 Obs_Infl*Obs_r40;
 Obs_Infl*Obs_xhr40;
-%Obs_Infl*Obs_GoY;
-%Obs_Infl*Obs_hours;
 
+Obs_r1*Obs_r1;
 Obs_r1*Obs_r40;
 Obs_r1*Obs_xhr40;
-%Obs_r1*Obs_GoY;
-%Obs_r1*Obs_hours;
 
+Obs_r40*Obs_r40;
 Obs_r40*Obs_xhr40;
-%Obs_r40*Obs_GoY;
-%Obs_r40*Obs_hours;
 
-%Obs_xhr40*Obs_GoY;
-%Obs_xhr40*Obs_hours;
+Obs_xhr40*Obs_xhr40;
 
-%Obs_GoY*Obs_hours;
+Obs_GoY*Obs_GoY;
+
+Obs_hours*Obs_hours;
 
-%first autocovariance
 Obs_Gr_C*Obs_Gr_C(-1);
 Obs_Gr_I*Obs_Gr_I(-1);
 Obs_Infl*Obs_Infl(-1);
@@ -513,6 +499,7 @@ Obs_GoY*Obs_GoY(-1);
 Obs_hours*Obs_hours(-1);
 end;
 
+
 %--------------------------------------------------------------------------
 % Create Data
 %--------------------------------------------------------------------------

tests/estimation/method_of_moments/AnScho/AnScho_GMM_estimParams0.mod

@@ -0,0 +1,5 @@
+% Test estimated_params block using initial values without bounds
+
+@#define estimParams = 0
+@#define MoM_Method = "GMM"
+@#include "AnScho_MoM_common.inc"

tests/estimation/method_of_moments/AnScho/AnScho_GMM_estimParams1.mod

@@ -0,0 +1,5 @@
+% Test estimated_params block using initial values and bounds
+
+@#define estimParams = 1
+@#define MoM_Method = "GMM"
+@#include "AnScho_MoM_common.inc"

tests/estimation/method_of_moments/AnScho/AnScho_GMM_estimParams2.mod

@@ -0,0 +1,5 @@
+% Test estimated_params block using prior information
+
+@#define estimParams = 2
+@#define MoM_Method = "GMM"
+@#include "AnScho_MoM_common.inc"

tests/estimation/method_of_moments/AnScho/AnScho_MoM.mod

@@ -1,286 +0,0 @@
-% DSGE model used in replication files of
-% An, Sungbae and Schorfheide, Frank, (2007), Bayesian Analysis of DSGE Models, Econometric Reviews, 26, issue 2-4, p. 113-172.
-% Adapted by Willi Mutschler (@wmutschl, willi@mutschler.eu)
-% =========================================================================
-% Copyright (C) 2020-2021 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 <https://www.gnu.org/licenses/>.
-% =========================================================================
-
-% Define testscenario
-@#define orderApp = 2
-@#define estimParams = 1
-
-% Note that we set the numerical optimization tolerance levels very large to speed up the testsuite
-@#define optimizer = 13
-
-var c p R g y z INFL INT YGR;
-varexo e_r e_g e_z;
-parameters tau nu kap cyst psi1 psi2 rhor rhog rhoz rrst pist gamst;
-
-varobs INT YGR INFL;
-
-tau   = 2;
-nu    = 0.1;
-kap   = 0.33;
-cyst  = 0.85;
-psi1  = 1.5;
-psi2  = 0.125;
-rhor  = 0.75;
-rhog  = 0.95;
-rhoz  = 0.9;
-rrst  = 1;
-pist  = 3.2;
-gamst = 0.55;
-
-model;
-#pist2 = exp(pist/400);
-#rrst2 = exp(rrst/400);
-#bet   = 1/rrst2;
-#phi   = tau*(1-nu)/nu/kap/pist2^2;
-#gst   = 1/cyst;
-#cst   = (1-nu)^(1/tau);
-#yst   = cst*gst;
-#dy    = y-y(-1);
-1 = exp(-tau*c(+1)+tau*c+R-z(+1)-p(+1));
-(1-nu)/nu/phi/(pist2^2)*(exp(tau*c)-1) = (exp(p)-1)*((1-1/2/nu)*exp(p)+1/2/nu) - bet*(exp(p(+1))-1)*exp(-tau*c(+1)+tau*c+y(+1)-y+p(+1));
-exp(c-y) = exp(-g) - phi*pist2^2*gst/2*(exp(p)-1)^2;
-R = rhor*R(-1) + (1-rhor)*psi1*p + (1-rhor)*psi2*(dy+z) + e_r/100;
-g = rhog*g(-1) + e_g/100;
-z = rhoz*z(-1) + e_z/100;
-YGR = gamst+100*(dy+z);
-INFL = pist+400*p;
-INT = pist+rrst+4*gamst+400*R;
-end;
-
-steady_state_model;
-  z = 0; p = 0; g = 0; r = 0; c = 0; y = 0;
-  YGR = gamst; INFL = pist; INT = pist + rrst + 4*gamst;
-end;
-
-shocks;
-  var e_r = 0.20^2;
-  var e_g = 0.80^2;
-  var e_z = 0.45^2;
-  corr e_r,e_g = 0.2;
-end;
-
-@#if estimParams == 0
-% Define only initial values without bounds
-estimated_params;
-  %tau,             1.50;
-  %kap,             0.15;
-  psi1,            1.20;
-  psi2,            0.50;
-  rhor,            0.50;
-  %rhog,            0.50;
-  %rhoz,            0.50;
-  %rrst,            1.20;
-  %pist,            3.00;
-  gamst,           0.75;
-  stderr e_r,      0.30;
-  stderr e_g,      0.30;
-  stderr e_z,      0.30;
-  corr e_r,e_g,    0.10;
-end;
-@#endif
-
-@#if estimParams == 1
-% Define initial values and bounds
-estimated_params;
-  %tau,             1.50,         1e-5,        10;
-  %kap,             0.15,         1e-5,        10;
-  psi1,            1.20,         1e-5,        10;
-  psi2,            0.50,         1e-5,        10;
-  rhor,            0.50,         1e-5,        0.99999;
-  %rhog,            0.50,         1e-5,        0.99999;
-  %rhoz,            0.50,         1e-5,        0.99999;
-  %rrst,            1.20,         1e-5,        10;
-  %pist,            3.00,         1e-5,        20;
-  gamst,           0.75,         -5,          5;
-  stderr e_r,      0.30,         1e-8,        5;
-  stderr e_g,      0.30,         1e-8,        5;
-  stderr e_z,      0.30,         1e-8,        5;
-  corr e_r,e_g,    0.10,         -1,          1;
-end;
-@#endif
-
-@#if estimParams == 2
-% Define prior distribution
-estimated_params;
-  %tau,             1.50,          1e-5,        10,          gamma_pdf,     2.00,       0.50;
-  %kap,             0.15,          1e-5,        10,          gamma_pdf,     0.33,       0.10;
-  psi1,            1.20,          1e-5,        10,          gamma_pdf,     1.50,       0.25;
-  psi2,            0.50,          1e-5,        10,          gamma_pdf,     0.125,      0.25;
-  rhor,            0.50,          1e-5,        0.99999,     beta_pdf,      0.50,       0.20;
-  %rhog,            0.50,          1e-5,        0.99999,     beta_pdf,      0.80,       0.10;
-  %rhoz,            0.50,          1e-5,        0.99999,     beta_pdf,      0.66,       0.15;
-  %rrst,            1.20,          1e-5,        10,          gamma_pdf,     0.50,       0.50;
-  %pist,            3.00,          1e-5,        20,          gamma_pdf,     7.00,       2.00;
-  gamst,           0.75,          -5,          5,           normal_pdf,    0.40,       0.20;
-  stderr e_r,      0.30,          1e-8,        5,           inv_gamma_pdf, 0.50,       0.26;
-  stderr e_g,      0.30,          1e-8,        5,           inv_gamma_pdf, 1.25,       0.65;
-  stderr e_z,      0.30,          1e-8,        5,           inv_gamma_pdf, 0.63,       0.33;
-  corr e_r,e_g,    0.10,          -1,          1,           uniform_pdf,       ,           , -1, 1;
-end;
-@#endif
-
-
-% Simulate data
-stoch_simul(order=@{orderApp},pruning,nodisplay,nomoments,periods=750,drop=500);
-save('AnScho_MoM_data_@{orderApp}.mat', options_.varobs{:} );
-pause(1);
-
-
-%--------------------------------------------------------------------------
-% Method of Moments Estimation
-%--------------------------------------------------------------------------
-matched_moments;
-YGR;
-INFL;
-INT;
-%second-order contemporenous product moments
-YGR*YGR;
-YGR*INFL;
-YGR*INT;
-INFL*INFL;
-INFL*INT;
-INT*INT;
-%second-order temporal product moments
-YGR*YGR(-1);
-INT*INT(-1);
-INFL*INFL(-1);
-end;
-
-% get indices in declaration order
-iYGR  = strmatch('YGR',  M_.endo_names,'exact');
-iINFL = strmatch('INFL', M_.endo_names,'exact');
-iINT  = strmatch('INT',  M_.endo_names,'exact');
-% first entry: number of variable in declaration order
-% second entry: lag
-% third entry: power
-
-matched_moments_ = {
-    %first-order product moments
-    [iYGR       ]  [0   ],  [1  ];
-    [iINFL      ]  [0   ],  [1  ];
-    [iINT       ]  [0   ],  [1  ];
-    %second-order contemporenous product moments
-    [iYGR  iYGR ]  [0  0],  [1 1];
-    [iYGR  iINFL]  [0  0],  [1 1];
-    [iYGR  iINT ]  [0  0],  [1 1];
-    [iINFL iINFL]  [0  0],  [1 1];
-    [iINFL iINT ]  [0  0],  [1 1];
-    [iINT  iINT ]  [0  0],  [1 1];
-    %second-order temporal product moments
-    [iYGR  iYGR ]  [0 -1],  [1 1];
-    %[iINT  iYGR ]  [0 -1],  [1 1];
-    %[iINFL iYGR ]  [0 -1],  [1 1];
-    %[iYGR  iINT ]  [0 -1],  [1 1];
-    [iINT  iINT ]  [0 -1],  [1 1];
-    %[iINFL iINT ]  [0 -1],  [1 1];
-    %[iYGR  iINFL]  [0 -1],  [1 1];
-    %[iINT  iINFL]  [0 -1],  [1 1];
-    [iINFL iINFL]  [0 -1],  [1 1];
-};
-
-if ~isequal(M_.matched_moments,matched_moments_)
-    error('Translation to matched_moments-block failed')
-end
-
-@#for mommethod in ["GMM", "SMM"]
-    method_of_moments(
-    % Necessery options
-          mom_method = @{mommethod}           % method of moments method; possible values: GMM|SMM
-        , datafile   = 'AnScho_MoM_data_@{orderApp}.mat' % name of filename with data
-
-    % Options for both GMM and SMM
-        % , bartlett_kernel_lag = 20          % bandwith in optimal weighting matrix
-        , order = @{orderApp}                 % order of Taylor approximation in perturbation
-        % , penalized_estimator               % include deviation from prior mean as additional moment restriction and use prior precision as weight
-        , pruning                             % use pruned state space system at higher-order
-        % , verbose                           % display and store intermediate estimation results
-        , weighting_matrix = ['optimal']      % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename. Size of cell determines stages in iterated estimation, e.g. two state with ['DIAGONAL','OPTIMAL']
-        %, weighting_matrix_scaling_factor=1  % scaling of weighting matrix in objective function
-        , se_tolx=1e-6                        % step size for numerical computation of standard errors
-
-    % Options for SMM
-        % , burnin=500                        % number of periods dropped at beginning of simulation
-        % , bounded_shock_support             % trim shocks in simulation to +- 2 stdev
-        % , seed = 24051986                   % seed used in simulations
-        % , simulation_multiple = 5           % multiple of the data length used for simulation
-
-    % Options for GMM
-        @#if mommethod == "GMM"
-        , analytic_standard_errors            % compute standard errors using analytical derivatives
-        @#endif
-
-    % General options
-        % , dirname = 'MM'                    % directory in which to store estimation output
-        % , graph_format = EPS                % specify the file format(s) for graphs saved to disk
-        % , nodisplay                         % do not display the graphs, but still save them to disk
-        % , nograph                           % do not create graphs (which implies that they are not saved to the disk nor displayed)
-        % , noprint                           % do not print stuff to console
-        % , plot_priors = 1                   % control plotting of priors
-        % , prior_trunc = 1e-10               % probability of extreme values of the prior density that is ignored when computing bounds for the parameters
-        % , TeX                               % print TeX tables and graphics
-
-    % Data and model options
-        % , first_obs = 501                   % number of first observation
-        % , logdata                           % if data is already in logs
-        , nobs = 250                          % number of observations
-        % , prefilter=0                       % demean each data series by its empirical mean and use centered moments
-
-        % , xls_sheet = data                  % name/number of sheet with data in Excel
-        % , xls_range = B2:D200               % range of data in Excel sheet
-
-    % Optimization options that can be set by the user in the mod file, otherwise default values are provided
-        % , huge_number=1e7                   % value for replacing the infinite bounds on parameters by finite numbers. Used by some optimizers for numerical reasons
-        , mode_compute = @{optimizer}         % specifies the optimizer for minimization of moments distance
-        , additional_optimizer_steps = [1]    % vector of additional mode-finders run after mode_compute
-        % optim: a list of NAME and VALUE pairs to set options for the optimization routines. Available options depend on mode_compute, some exemplary common options:
-        , optim = ('TolFun'      , 1e-6       % termination tolerance on the function value, a positive scalar
-                  ,'TolX'        , 1e-6       % termination tolerance on x, a positive scalar
-                  ,'MaxIter'     , 3000       % maximum number of iterations allowed, a positive integer
-                  ,'MaxFunEvals' , 1D6        % maximum number of function evaluations allowed, a positive integer
-        %           ,'UseParallel' , 1        % when true (and supported by optimizer) solver estimates gradients in parallel (using Matlab/Octave's parallel toolbox)
-        %           ,'Jacobian'    , 'off'    % when 'off' gradient-based solvers approximate Jacobian using finite differences; for GMM we can also pass the analytical Jacobian to gradient-based solvers by setting this 'on'
-                  )                         
-        , silent_optimizer                    % run minimization of moments distance silently without displaying results or saving files in between
-
-    % Numerical algorithms options
-        % , aim_solver                             % Use AIM algorithm to compute perturbation approximation
-        % , k_order_solver                         % use k_order_solver in higher order perturbation approximations
-        % , dr=default                             % method used to compute the decision rule; possible values are DEFAULT, CYCLE_REDUCTION, LOGARITHMIC_REDUCTION
-        % , dr_cycle_reduction_tol = 1e-7          % convergence criterion used in the cycle reduction algorithm
-        % , dr_logarithmic_reduction_tol = 1e-12   % convergence criterion used in the logarithmic reduction algorithm
-        % , dr_logarithmic_reduction_maxiter = 100 % maximum number of iterations used in the logarithmic reduction algorithm
-        % , lyapunov = DEFAULT                     % algorithm used to solve lyapunov equations; possible values are DEFAULT, FIXED_POINT, DOUBLING, SQUARE_ROOT_SOLVER
-        % , lyapunov_complex_threshold = 1e-15     % complex block threshold for the upper triangular matrix in symmetric Lyapunov equation solver
-        % , lyapunov_fixed_point_tol = 1e-10       % convergence criterion used in the fixed point Lyapunov solver
-        % , lyapunov_doubling_tol = 1e-16          % convergence criterion used in the doubling algorithm
-        % , sylvester = default                    % algorithm to solve Sylvester equation; possible values are DEFAULT, FIXED_POINT
-        % , sylvester_fixed_point_tol = 1e-12      % convergence criterion used in the fixed point Sylvester solver
-        % , qz_criterium = 0.999999                % value used to split stable from unstable eigenvalues in reordering the Generalized Schur decomposition used for solving first order problems
-        % , qz_zero_threshold = 1e-6               % value used to test if a generalized eigenvalue is 0/0 in the generalized Schur decomposition
-        % , schur_vec_tol=1e-11                    % tolerance level used to find nonstationary variables in Schur decomposition of the transition matrix
-        % , mode_check                             % plot the target function for values around the computed minimum for each estimated parameter in turn
-        % , mode_check_neighbourhood_size = 5      % width of the window (expressed in percentage deviation) around the computed minimum to be displayed on the diagnostic plots
-        % , mode_check_symmetric_plots=1           % ensure that the check plots are symmetric around the minimum
-        % , mode_check_number_of_points = 20       % number of points around the minimum where the target function is evaluated (for each parameter)
-    );
-@#endfor
-

tests/estimation/method_of_moments/AnScho/AnScho_MoM_common.inc

@@ -0,0 +1,279 @@
+% DSGE model used in replication files of
+% An, Sungbae and Schorfheide, Frank, (2007), Bayesian Analysis of DSGE Models, Econometric Reviews, 26, issue 2-4, p. 113-172.
+% Adapted by Willi Mutschler (@wmutschl, willi@mutschler.eu)
+% =========================================================================
+% Copyright (C) 2020-2021 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 <https://www.gnu.org/licenses/>.
+% =========================================================================
+
+var c p R g y z INFL INT YGR;
+varexo e_r e_g e_z;
+parameters tau nu kap cyst psi1 psi2 rhor rhog rhoz rrst pist gamst;
+
+varobs INT YGR INFL;
+
+tau   = 2;
+nu    = 0.1;
+kap   = 0.33;
+cyst  = 0.85;
+psi1  = 1.5;
+psi2  = 0.125;
+rhor  = 0.75;
+rhog  = 0.95;
+rhoz  = 0.9;
+rrst  = 1;
+pist  = 3.2;
+gamst = 0.55;
+
+model;
+#pist2 = exp(pist/400);
+#rrst2 = exp(rrst/400);
+#bet   = 1/rrst2;
+#phi   = tau*(1-nu)/nu/kap/pist2^2;
+#gst   = 1/cyst;
+#cst   = (1-nu)^(1/tau);
+#yst   = cst*gst;
+#dy    = y-y(-1);
+1 = exp(-tau*c(+1)+tau*c+R-z(+1)-p(+1));
+(1-nu)/nu/phi/(pist2^2)*(exp(tau*c)-1) = (exp(p)-1)*((1-1/2/nu)*exp(p)+1/2/nu) - bet*(exp(p(+1))-1)*exp(-tau*c(+1)+tau*c+y(+1)-y+p(+1));
+exp(c-y) = exp(-g) - phi*pist2^2*gst/2*(exp(p)-1)^2;
+R = rhor*R(-1) + (1-rhor)*psi1*p + (1-rhor)*psi2*(dy+z) + e_r/100;
+g = rhog*g(-1) + e_g/100;
+z = rhoz*z(-1) + e_z/100;
+YGR = gamst+100*(dy+z);
+INFL = pist+400*p;
+INT = pist+rrst+4*gamst+400*R;
+end;
+
+steady_state_model;
+  z = 0; p = 0; g = 0; R = 0; c = 0; y = 0;
+  YGR = gamst; INFL = pist; INT = pist + rrst + 4*gamst;
+end;
+
+shocks;
+  var e_r = 0.20^2;
+  var e_g = 0.80^2;
+  var e_z = 0.45^2;
+  corr e_r,e_g = 0.2;
+end;
+
+@#if estimParams == 0
+% Define only initial values without bounds
+estimated_params;
+  %tau,             1.50;
+  %kap,             0.15;
+  psi1,            1.20;
+  psi2,            0.50;
+  rhor,            0.50;
+  %rhog,            0.50;
+  %rhoz,            0.50;
+  %rrst,            1.20;
+  %pist,            3.00;
+  gamst,           0.75;
+  stderr e_r,      0.30;
+  stderr e_g,      0.30;
+  stderr e_z,      0.30;
+  corr e_r,e_g,    0.10;
+end;
+@#endif
+
+@#if estimParams == 1
+% Define initial values and bounds
+estimated_params;
+  %tau,             1.50,         1e-5,        10;
+  %kap,             0.15,         1e-5,        10;
+  psi1,            1.20,         1e-5,        10;
+  psi2,            0.50,         1e-5,        10;
+  rhor,            0.50,         1e-5,        0.99999;
+  %rhog,            0.50,         1e-5,        0.99999;
+  %rhoz,            0.50,         1e-5,        0.99999;
+  %rrst,            1.20,         1e-5,        10;
+  %pist,            3.00,         1e-5,        20;
+  gamst,           0.75,         -5,          5;
+  stderr e_r,      0.30,         1e-8,        5;
+  stderr e_g,      0.30,         1e-8,        5;
+  stderr e_z,      0.30,         1e-8,        5;
+  corr e_r,e_g,    0.10,         -1,          1;
+end;
+@#endif
+
+@#if estimParams == 2
+% Define prior distribution
+estimated_params;
+  %tau,             1.50,          1e-5,        10,          gamma_pdf,     2.00,       0.50;
+  %kap,             0.15,          1e-5,        10,          gamma_pdf,     0.33,       0.10;
+  psi1,            1.20,          1e-5,        10,          gamma_pdf,     1.50,       0.25;
+  psi2,            0.50,          1e-5,        10,          gamma_pdf,     0.125,      0.25;
+  rhor,            0.50,          1e-5,        0.99999,     beta_pdf,      0.50,       0.20;
+  %rhog,            0.50,          1e-5,        0.99999,     beta_pdf,      0.80,       0.10;
+  %rhoz,            0.50,          1e-5,        0.99999,     beta_pdf,      0.66,       0.15;
+  %rrst,            1.20,          1e-5,        10,          gamma_pdf,     0.50,       0.50;
+  %pist,            3.00,          1e-5,        20,          gamma_pdf,     7.00,       2.00;
+  gamst,           0.75,          -5,          5,           normal_pdf,    0.40,       0.20;
+  stderr e_r,      0.30,          1e-8,        5,           inv_gamma_pdf, 0.50,       0.26;
+  stderr e_g,      0.30,          1e-8,        5,           inv_gamma_pdf, 1.25,       0.65;
+  stderr e_z,      0.30,          1e-8,        5,           inv_gamma_pdf, 0.63,       0.33;
+  corr e_r,e_g,    0.10,          -1,          1,           uniform_pdf,       ,           , -1, 1;
+end;
+@#endif
+
+
+% % Simulate data
+% stoch_simul(order=2,pruning,nodisplay,nomoments,periods=750,drop=500);
+% save('AnScho_MoM_data_2.mat', options_.varobs{:} );
+% pause(1);
+
+%--------------------------------------------------------------------------
+% Method of Moments Estimation
+%--------------------------------------------------------------------------
+matched_moments;
+
+%first-order product moments
+YGR;
+YGR(-1);   %redundant
+YGR(0)^1;  %redundant
+YGR^1;     %redundant
+YGR^1*INFL^0*INT^0; %redundant
+INFL(+2)^1;
+INT(-2);
+INT(2);    %redundant
+
+%second-order contemporenous product moments
+YGR^2;
+YGR(-1)^2; %redundant
+YGR*YGR;   %redundant
+YGR(1)*YGR(1);  %redundant
+INFL*YGR;
+YGR*INFL;  %redundant
+YGR(2)*INT(2);
+INT(2)*YGR(2); %redundant
+INT(2)^1*YGR(2)^1; %redundant
+YGR(2)^1*INT(2)^1; %redundant
+INFL(-2)^1*INFL(-2)^1;
+INFL(-1)*INT(-1);
+INT(0)^1*INT^1;
+INT(0)^1*INT^1*YGR(2)^0*INT(-2)^0*INFL^0; %redundant
+
+%second-order temporal product moments
+YGR*YGR(-1);
+YGR(3)*YGR(2); %redundant
+YGR(2)*YGR(3); %redundant
+YGR(-2)*YGR(-1); %redundant
+INT(3)^1*INT(-2)^1;
+YGR(5)^0*INFL*INFL(2)^1;
+YGR(5)^1*INFL(-3)^1;
+INFL(-3)^1*YGR(5)^1; %redundant
+INFL(3)*INT(2);
+
+@#ifdef Extended_Matched_Moments_Checks
+YGR^3;
+YGR(-3)^2*YGR(-3); %redundant
+YGR(-3)*YGR(-3)^2; %redundant
+YGR(0)^1*YGR(0)*YGR; %redundant
+INT(-2)^2*INFL(-1)^4;
+INFL(1)^4*INT(0)^2; %redundant
+INT(0)^2*INFL(1)^4; %redundant
+YGR^2*INFL(-3)^4*INT(5)^6;
+YGR^2*INT(5)^6*INFL(-3)^4;%redundant
+INT(5)^6*YGR^2*INFL(-3)^4;%redundant
+@#endif
+end;
+
+
+method_of_moments(
+% Necessery options
+      mom_method = @{MoM_Method}           % method of moments method; possible values: GMM|SMM
+    , datafile   = 'AnScho_MoM_data_2.mat' % name of filename with data
+
+% Options for both GMM and SMM
+    % , bartlett_kernel_lag = 20          % bandwith in optimal weighting matrix
+    , order = 2                           % order of Taylor approximation in perturbation
+    % , penalized_estimator               % include deviation from prior mean as additional moment restriction and use prior precision as weight
+    , pruning                             % use pruned state space system at higher-order
+    , verbose                           % display and store intermediate estimation results
+    , weighting_matrix = ['optimal']      % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename. Size of cell determines stages in iterated estimation, e.g. two state with ['DIAGONAL','OPTIMAL']
+    %, weighting_matrix_scaling_factor=1  % scaling of weighting matrix in objective function
+    , se_tolx=1e-6                        % step size for numerical computation of standard errors
+
+% Options for SMM
+        % , burnin=500                        % number of periods dropped at beginning of simulation
+        , bounded_shock_support             % trim shocks in simulation to +- 2 stdev
+        % , seed = 24051986                   % seed used in simulations
+        % , simulation_multiple = 5           % multiple of the data length used for simulation
+
+% Options for GMM
+@#if MoM_Method == "GMM"
+    , analytic_standard_errors            % compute standard errors using analytical derivatives
+@#endif
+% General options
+    % , dirname = 'MM'                    % directory in which to store estimation output
+    % , graph_format = EPS                % specify the file format(s) for graphs saved to disk
+    % , nodisplay                         % do not display the graphs, but still save them to disk
+    % , nograph                           % do not create graphs (which implies that they are not saved to the disk nor displayed)
+    % , noprint                           % do not print stuff to console
+    % , plot_priors = 1                   % control plotting of priors
+    % , prior_trunc = 1e-10               % probability of extreme values of the prior density that is ignored when computing bounds for the parameters
+    % , TeX                               % print TeX tables and graphics
+
+% Data and model options
+    % , first_obs = 501                   % number of first observation
+    % , logdata                           % if data is already in logs
+    , nobs = 250                          % number of observations
+    % , prefilter=0                       % demean each data series by its empirical mean and use centered moments
+
+    % , xls_sheet = data                  % name/number of sheet with data in Excel
+    % , xls_range = B2:D200               % range of data in Excel sheet
+
+% Optimization options that can be set by the user in the mod file, otherwise default values are provided
+    % , huge_number=1e7                   % value for replacing the infinite bounds on parameters by finite numbers. Used by some optimizers for numerical reasons
+    @#ifdef NoEstim
+    , mode_compute = 0
+    @#else
+    , mode_compute = 13                   % specifies the optimizer for minimization of moments distance
+    , additional_optimizer_steps = [1]    % vector of additional mode-finders run after mode_compute
+    @#endif
+    % optim: a list of NAME and VALUE pairs to set options for the optimization routines. Available options depend on mode_compute, some exemplary common options:
+    , optim = ('TolFun'      , 1e-6       % termination tolerance on the function value, a positive scalar
+              ,'TolX'        , 1e-6       % termination tolerance on x, a positive scalar
+              ,'MaxIter'     , 3000       % maximum number of iterations allowed, a positive integer
+              ,'MaxFunEvals' , 1D6        % maximum number of function evaluations allowed, a positive integer
+    %           ,'UseParallel' , 1        % when true (and supported by optimizer) solver estimates gradients in parallel (using Matlab/Octave's parallel toolbox)
+    %           ,'Jacobian'    , 'off'    % when 'off' gradient-based solvers approximate Jacobian using finite differences; for GMM we can also pass the analytical Jacobian to gradient-based solvers by setting this 'on'
+              )
+    %, silent_optimizer                    % run minimization of moments distance silently without displaying results or saving files in between
+
+% Numerical algorithms options
+    % , aim_solver                             % Use AIM algorithm to compute perturbation approximation
+    % , k_order_solver                         % use k_order_solver in higher order perturbation approximations
+    % , dr=default                             % method used to compute the decision rule; possible values are DEFAULT, CYCLE_REDUCTION, LOGARITHMIC_REDUCTION
+    % , dr_cycle_reduction_tol = 1e-7          % convergence criterion used in the cycle reduction algorithm
+    % , dr_logarithmic_reduction_tol = 1e-12   % convergence criterion used in the logarithmic reduction algorithm
+    % , dr_logarithmic_reduction_maxiter = 100 % maximum number of iterations used in the logarithmic reduction algorithm
+    % , lyapunov = DEFAULT                     % algorithm used to solve lyapunov equations; possible values are DEFAULT, FIXED_POINT, DOUBLING, SQUARE_ROOT_SOLVER
+    % , lyapunov_complex_threshold = 1e-15     % complex block threshold for the upper triangular matrix in symmetric Lyapunov equation solver
+    % , lyapunov_fixed_point_tol = 1e-10       % convergence criterion used in the fixed point Lyapunov solver
+    % , lyapunov_doubling_tol = 1e-16          % convergence criterion used in the doubling algorithm
+    % , sylvester = default                    % algorithm to solve Sylvester equation; possible values are DEFAULT, FIXED_POINT
+    % , sylvester_fixed_point_tol = 1e-12      % convergence criterion used in the fixed point Sylvester solver
+    % , qz_criterium = 0.999999                % value used to split stable from unstable eigenvalues in reordering the Generalized Schur decomposition used for solving first order problems
+    % , qz_zero_threshold = 1e-6               % value used to test if a generalized eigenvalue is 0/0 in the generalized Schur decomposition
+    % , schur_vec_tol=1e-11                    % tolerance level used to find nonstationary variables in Schur decomposition of the transition matrix
+    , mode_check                             % plot the target function for values around the computed minimum for each estimated parameter in turn
+    % , mode_check_neighbourhood_size = 5      % width of the window (expressed in percentage deviation) around the computed minimum to be displayed on the diagnostic plots
+    % , mode_check_symmetric_plots=1           % ensure that the check plots are symmetric around the minimum
+    % , mode_check_number_of_points = 20       % number of points around the minimum where the target function is evaluated (for each parameter)
+);

tests/estimation/method_of_moments/AnScho/AnScho_SMM_estimParams0.mod

@@ -0,0 +1,5 @@
+% Test estimated_params block using initial values without bounds
+
+@#define estimParams = 0
+@#define MoM_Method = "SMM"
+@#include "AnScho_MoM_common.inc"

tests/estimation/method_of_moments/AnScho/AnScho_SMM_estimParams1.mod

@@ -0,0 +1,5 @@
+% Test estimated_params block using initial values and bounds
+
+@#define estimParams = 1
+@#define MoM_Method = "SMM"
+@#include "AnScho_MoM_common.inc"

tests/estimation/method_of_moments/AnScho/AnScho_SMM_estimParams2.mod

@@ -0,0 +1,5 @@
+% Test estimated_params block using prior information
+
+@#define estimParams = 2
+@#define MoM_Method = "SMM"
+@#include "AnScho_MoM_common.inc"

tests/estimation/method_of_moments/AnScho/AnScho_matched_moments.mod

@@ -0,0 +1,134 @@
+% Test translation from matched_moments block to M_.matched_moments
+
+@#define estimParams = 0
+@#define MoM_Method = "SMM"
+@#define NoEstim = 1
+@#define Extended_Matched_Moments_Checks
+@#include "AnScho_MoM_common.inc"
+
+% get indices in declaration order
+iYGR  = strmatch('YGR',  M_.endo_names,'exact');
+iINFL = strmatch('INFL', M_.endo_names,'exact');
+iINT  = strmatch('INT',  M_.endo_names,'exact');
+
+% M_.matched_moments has the following structure:
+%   - first entry: index number of variable in declaration order
+%   - second entry: lead or lag
+%   - third entry: power
+
+matched_moments_orig = {
+
+    %first-order product moments
+    [iYGR       ]  [0 ],  [1];
+    [iYGR       ]  [-1],  [1];
+    [iYGR       ]  [0 ],  [1];
+    [iYGR       ]  [0 ],  [1];
+    [iYGR       ]  [0 ],  [1];
+    [iINFL      ]  [2 ],  [1];
+    [iINT       ]  [-2],  [1];
+    [iINT       ]  [2 ],  [1];
+
+    %second-order contemporenous product moments
+    [iYGR       ]  [0     ],  [2  ];
+    [iYGR       ]  [-1    ],  [2  ];
+    [iYGR  iYGR ]  [0    0],  [1 1];
+    [iYGR  iYGR ]  [1    1],  [1 1];
+    [iYGR  iINFL]  [0    0],  [1 1];
+    [iYGR  iINFL]  [0    0],  [1 1];
+    [iINT  iYGR ]  [2    2],  [1 1];
+    [iINT  iYGR ]  [2    2],  [1 1];
+    [iINT  iYGR ]  [2    2],  [1 1];
+    [iINT  iYGR ]  [2    2],  [1 1];
+    [iINFL iINFL]  [-2  -2],  [1 1];
+    [iINFL iINT ]  [-1  -1],  [1 1];
+    [iINT  iINT ]  [0    0],  [1 1];
+    [iINT  iINT ]  [0    0],  [1 1];
+
+    %second-order temporal product moments
+    [iYGR  iYGR ]  [0  -1],  [1 1];
+    [iYGR  iYGR ]  [2   3],  [1 1];
+    [iYGR  iYGR ]  [2   3],  [1 1];
+    [iYGR  iYGR ]  [-1 -2],  [1 1];
+    [iINT  iINT ]  [-2  3],  [1 1];
+    [iINFL iINFL]  [0   2],  [1 1];
+    [iYGR  iINFL]  [5  -3],  [1 1];
+    [iYGR  iINFL]  [5  -3],  [1 1];
+    [iINT  iINFL]  [2   3],  [1 1];
+
+    [iYGR            ]  [0      ],  [3     ];
+    [iYGR  iYGR      ]  [-3  -3 ],  [1  2  ];
+    [iYGR  iYGR      ]  [-3  -3 ],  [1  2  ];
+    [iYGR  iYGR  iYGR]  [0  0  0],  [1  1  1];
+    [iINT  iINFL     ]  [-2  -1 ],  [2  4  ];
+    [iINFL iINT      ]  [1   0 ],   [4  2  ];
+    [iINFL iINT      ]  [1   0 ],   [4  2  ];
+    [iYGR  iINFL iINT]  [0  -3  5],   [2  4  6];
+    [iINFL iYGR  iINT]  [-3  0  5],   [4  2  6];
+    [iINFL iYGR  iINT]  [-3  0  5],   [4  2  6];
+};
+
+% Removed duplicate moment conditions
+matched_moments_no_duplicate= {
+
+    %first-order product moments
+    [iYGR       ]  [0 ],  [1];
+%    [iYGR       ]  [-1],  [1];
+%    [iYGR       ]  [0 ],  [1];
+%    [iYGR       ]  [0 ],  [1];
+%    [iYGR       ]  [0 ],  [1];
+    [iINFL      ]  [2 ],  [1];
+    [iINT       ]  [-2],  [1];
+%    [iINT       ]  [2 ],  [1];
+
+    %second-order contemporenous product moments
+    [iYGR       ]  [0     ],  [2  ];
+%    [iYGR       ]  [-1    ],  [2  ];
+%    [iYGR  iYGR ]  [0    0],  [1 1];
+%    [iYGR  iYGR ]  [1    1],  [1 1];
+    [iYGR  iINFL]  [0    0],  [1 1];
+%    [iYGR  iINFL]  [0    0],  [1 1];
+    [iINT  iYGR ]  [2    2],  [1 1];
+%    [iINT  iYGR ]  [2    2],  [1 1];
+%    [iINT  iYGR ]  [2    2],  [1 1];
+%    [iINT  iYGR ]  [2    2],  [1 1];
+    [iINFL iINFL]  [-2  -2],  [1 1];
+    [iINFL iINT ]  [-1  -1],  [1 1];
+    [iINT  iINT ]  [0    0],  [1 1];
+%    [iINT  iINT ]  [0    0],  [1 1];
+
+    %second-order temporal product moments
+    [iYGR  iYGR ]  [0  -1],  [1 1];
+%    [iYGR  iYGR ]  [2   3],  [1 1];
+%    [iYGR  iYGR ]  [2   3],  [1 1];
+%    [iYGR  iYGR ]  [-1 -2],  [1 1];
+    [iINT  iINT ]  [-2  3],  [1 1];
+    [iINFL iINFL]  [0   2],  [1 1];
+    [iYGR  iINFL]  [5  -3],  [1 1];
+%    [iYGR  iINFL]  [5  -3],  [1 1];
+    [iINT  iINFL]  [2   3],  [1 1];
+
+    [iYGR            ]  [0      ],  [3     ];
+%    [iYGR  iYGR      ]  [-3  -3 ],  [1  2  ];
+%    [iYGR  iYGR      ]  [-3  -3 ],  [1  2  ];
+%    [iYGR  iYGR  iYGR]  [0  0  0],  [1  1  1];
+    [iINT  iINFL     ]  [-2  -1 ],  [2  4  ];
+%    [iINFL iINT      ]  [1   0 ],   [4  2  ];
+%    [iINFL iINT      ]  [1   0 ],   [4  2  ];
+    [iYGR  iINFL iINT]  [0  -3  5],   [2  4  6];
+%    [iINFL iYGR  iINT]  [-3  0  5],   [4  2  6];
+%    [iINFL iYGR  iINT]  [-3  0  5],   [4  2  6];
+};
+
+
+
+if ~isequal(M_.matched_moments_orig,matched_moments_orig)
+    error('Translation to matched_moments-block failed!')
+else
+    fprintf('Translation to matched_moments-block successful!\n\n')
+end
+
+if ~isequal(M_.matched_moments,matched_moments_no_duplicate)
+    error('Removal of duplicate moment conditions failed!')
+else
+    fprintf('Removal of duplicate moment conditions was successful!\n\n')
+end

tests/estimation/system_prior_restriction/Gali_2015.mod

@@ -0,0 +1,184 @@
+/*
+ * This file shows how to use "system prior"-type prior restrictions as in
+ * Michal Andrle/Miroslav Plašil (2018): "Econometrics with system priors",
+ * Economics Letters, 172, pp. 134-137 during estimation based on
+ * the baseline New Keynesian model of Jordi Galí (2015): Monetary Policy, Inflation,
+ * and the Business Cycle, Princeton University Press, Second Edition, Chapter 3
+ *
+ * THIS MOD-FILE REQUIRES DYNARE 4.5 OR HIGHER
+ *
+ * Notes:
+ *  - The estimation will automatically take the Gali_2015_prior_restrictions.m into
+ *      account, which has the required name and format
+ *  - Estimation is based on simulated data
+ *  - The file also shows how to use a prior/posterior-function
+ *
+ * This implementation was written by Johannes Pfeifer. In case you spot mistakes,
+ * email me at jpfeifer@gmx.de
+ *
+ * Please note that the following copyright notice only applies to this Dynare
+ */
+
+/*
+ * Copyright (C) 2021 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/>.
+ */
+
+var pi          ${\pi}$                 (long_name='inflation')
+    y_gap       ${\tilde y}$            (long_name='output gap')
+    y_nat       ${y^{nat}}$             (long_name='natural output')      //(in contrast to the textbook defined in deviation from steady state)
+    y           ${y}$                   (long_name='output')
+    yhat        ${\hat y}$              (long_name='output deviation from steady state')
+    r_nat       ${r^{nat}}$             (long_name='natural interest rate')
+    i           ${i}$                   (long_name='nominal interrst rate')
+    n           ${n}$                   (long_name='hours worked')
+    nu          ${\nu}$                 (long_name='AR(1) monetary policy shock process')
+    a           ${a}$                   (long_name='AR(1) technology shock process')
+    z           ${z}$                   (long_name='AR(1) preference shock process')
+    p           ${p}$                   (long_name='price level')
+    w           ${w}$                   (long_name='nominal wage')
+    c           ${c}$                   (long_name='consumption')
+;
+
+varexo  eps_a       ${\varepsilon_a}$       (long_name='technology shock')
+        eps_nu  ${\varepsilon_\nu}$     (long_name='monetary policy shock')
+        eps_z       ${\varepsilon_z}$   (long_name='preference shock innovation')
+       ;
+
+parameters alppha       ${\alpha}$     (long_name='capital share')
+    betta               ${\beta}$       (long_name='discount factor')
+    rho_a               ${\rho_a}$      (long_name='autocorrelation technology shock')
+    rho_nu          ${\rho_{\nu}}$  (long_name='autocorrelation monetary policy shock')
+    rho_z               ${\rho_{z}}$    (long_name='autocorrelation monetary demand shock')
+    siggma              ${\sigma}$      (long_name='inverse EIS')
+    varphi              ${\varphi}$     (long_name='inverse Frisch elasticity')
+    phi_pi              ${\phi_{\pi}}$  (long_name='inflation feedback Taylor Rule')
+    phi_y               ${\phi_{y}}$    (long_name='output feedback Taylor Rule')
+    eta                 ${\eta}$        (long_name='semi-elasticity of money demand')
+    epsilon             ${\epsilon}$    (long_name='demand elasticity')
+    theta               ${\theta}$      (long_name='Calvo parameter')
+    ;
+%----------------------------------------------------------------
+% Parametrization, p. 67  and p. 113-115
+%----------------------------------------------------------------
+siggma = 1;
+varphi = 5;
+phi_pi = 1.5;
+phi_y  = 0.125;
+theta = 3/4;
+rho_nu =0.5;
+rho_z  = 0.5;
+rho_a  = 0.9;
+betta  = 0.99;
+eta  = 3.77; %footnote 11, p. 115
+alppha = 1/4;
+epsilon = 9;
+
+%----------------------------------------------------------------
+% First Order Conditions
+%----------------------------------------------------------------
+
+model(linear);
+//Composite parameters
+#Omega=(1-alppha)/(1-alppha+alppha*epsilon);        %defined on page 60
+#psi_n_ya=(1+varphi)/(siggma*(1-alppha)+varphi+alppha);   %defined on page 62
+#lambda=(1-theta)*(1-betta*theta)/theta*Omega;      %defined on page 61
+#kappa=lambda*(siggma+(varphi+alppha)/(1-alppha));     %defined on page 63
+[name='New Keynesian Phillips Curve eq. (22)']
+pi=betta*pi(+1)+kappa*y_gap;
+[name='Dynamic IS Curve eq. (23)']
+y_gap=-1/siggma*(i-pi(+1)-r_nat)+y_gap(+1);
+[name='Interest Rate Rule eq. (26)']
+i=phi_pi*pi+phi_y*yhat+nu;
+[name='Definition natural rate of interest eq. (24)']
+r_nat=-siggma*psi_n_ya*(1-rho_a)*a+(1-rho_z)*z;
+[name='Definition natural output, eq. (20)']
+y_nat=psi_n_ya*a;
+[name='Definition output gap']
+y_gap=y-y_nat;
+[name='Monetary policy shock']
+nu=rho_nu*nu(-1)+eps_nu;
+[name='TFP shock']
+a=rho_a*a(-1)+eps_a;
+[name='Production function (eq. 14)']
+y=a+(1-alppha)*n;
+[name='Preference shock, p. 54']
+z=rho_z*z(-1) - eps_z;
+[name='Output deviation from steady state']
+yhat=y-steady_state(y);
+[name='Definition price level']
+pi=p-p(-1);
+[name='resource constraint, eq. (12)']
+y=c;
+[name='FOC labor, eq. (2)']
+w-p=siggma*c+varphi*n;
+end;
+
+
+shocks;
+    var eps_nu = 0.0025^2; //1 standard deviation shock of 25 basis points, i.e. 1 percentage point annualized
+end;
+
+% simulate data
+stoch_simul(periods=100,drop=0,irf=0) yhat;
+% save data
+datatomfile('sim_data',{'yhat'});
+
+estimated_params;
+theta,0.75,beta_pdf,0.5,0.1;
+betta, beta_pdf, 0.993, 0.002;
+alppha, beta_pdf, 0.25, 0.02;
+end;
+
+varobs yhat;
+
+% Run prior function to get prior slope of the PC based on independent priors
+hh=figure('Name','Slope of the Phillips Curve');
+prior_function(function='Gali_2015_PC_slope');
+PC_slope_vec=cell2mat(oo_.prior_function_results(:,1));
+
+optimal_bandwidth = mh_optimal_bandwidth(PC_slope_vec,length(PC_slope_vec),0,'gaussian');
+[density(:,1),density(:,2)] = kernel_density_estimate(PC_slope_vec,512,length(PC_slope_vec),optimal_bandwidth,'gaussian');
+figure(hh)
+subplot(3,1,1)
+plot(density(:,1),density(:,2));
+title('Prior')
+
+% Run estimation with 1 observation to show effect of _prior_restriction .m
+% on independent prior
+estimation(datafile='sim_data',mode_compute=5,mh_replic=2001,mh_nblocks=1,diffuse_filter,nobs=1,mh_jscale=0.8);
+posterior_function(function='Gali_2015_PC_slope');
+PC_slope_vec=cell2mat(oo_.posterior_function_results(:,1));
+optimal_bandwidth = mh_optimal_bandwidth(PC_slope_vec,length(PC_slope_vec),0,'gaussian');
+[density(:,1),density(:,2)] = kernel_density_estimate(PC_slope_vec,512,length(PC_slope_vec),optimal_bandwidth,'gaussian');
+figure(hh)
+subplot(3,1,2)
+plot(density(:,1),density(:,2));
+title('Updated Prior')
+
+
+% Run estimation with full observations
+estimation(datafile='sim_data',mode_compute=5,mh_replic=2001,mh_nblocks=1,diffuse_filter,nobs=100,mh_jscale=0.8);
+
+posterior_function(function='Gali_2015_PC_slope');
+PC_slope_vec=cell2mat(oo_.posterior_function_results(:,1));
+optimal_bandwidth = mh_optimal_bandwidth(PC_slope_vec,length(PC_slope_vec),0,'gaussian');
+[density(:,1),density(:,2)] = kernel_density_estimate(PC_slope_vec,512,length(PC_slope_vec),optimal_bandwidth,'gaussian');
+figure(hh)
+subplot(3,1,3)
+plot(density(:,1),density(:,2));
+title('Posterior')

tests/estimation/system_prior_restriction/Gali_2015_PC_slope.m

@@ -0,0 +1,56 @@
+function output_cell =PC_slope(xparam1,M_,options_,oo_,estim_params_,bayestopt_,dataset_,dataset_info)
+% output_cell =PC_slope(xparam1,M_,options_,oo_,estim_params_,bayestopt_,dataset_,dataset_info);
+% This is an example file computing statistics on the prior/posterior draws. The
+% function allows read-only access to all Dynare structures. However, those
+% structures are local to this function.  Changing them will not affect
+% other Dynare functions and you cannot use them to pass results to other
+% Dynare functions.
+% The function takes one and only one output argument: an 1 by n cell.
+% Using functions like cell2mat, the contents of the cell can be easily
+% transformed back to matrices. See the fs2000_posterior_function.mod for
+% an example
+
+% INPUTS
+%   xparam1                      Current parameter draw
+%   M_           [structure]     Matlab's structure describing the Model (initialized by dynare, see @ref{M_}).
+%   options_     [structure]     Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
+%   oo_          [structure]     Matlab's structure gathering the results (initialized by dynare, see @ref{oo_}).
+%   estim_params_[structure]     Matlab's structure describing the estimated_parameters (initialized by dynare, see @ref{estim_params_}).
+%   bayestopt_   [structure]     Matlab's structure describing the parameter options (initialized by dynare, see @ref{bayestopt_}).
+%   dataset_     [structure]     Matlab's structure storing the dataset
+%   dataset_info [structure]     Matlab's structure storing the information about the dataset
+
+% Output
+%   output_cell  [1 by n cell]   1 by n Matlab cell allowing to store any
+%                                desired computation or result (strings, matrices, structures, etc.)
+
+% Copyright (C) 2021 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/>.
+
+
+%% store the slope based on the parameter draw
+NumberOfParameters = M_.param_nbr;
+for ii = 1:NumberOfParameters
+  paramname = deblank(M_.param_names{ii,:});
+  eval([ paramname ' = M_.params(' int2str(ii) ');']);
+end
+
+Omega=(1-alppha)/(1-alppha+alppha*epsilon);
+lambda=(1-theta)*(1-betta*theta)/theta*Omega;      %defined on page 61
+output_cell{1,1}=lambda*(siggma+(varphi+alppha)/(1-alppha));     %defined on page 63
+
+end

tests/estimation/system_prior_restriction/Gali_2015_prior_restrictions.m

@@ -0,0 +1,40 @@
+function log_prior_val=Gali_2015_prior_restrictions(M_, oo_, options_, dataset_, dataset_info);
+% function prior_val=Gali_2015_prior_restrictions(M_, oo_, options_, dataset_, dataset_info);
+% Example of a _prior_restrictions-file automatically called during
+% estimation
+% It imposes a prior of the slope of the New Keynesian Phillips Curve of
+% 0.03. As the slope is a composite of other parameters with independent
+% priors, a separate function is required to do this.
+
+% Copyright (C) 2021 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 <https://www.gnu.org/licenses/>.
+
+
+% read out parameters to access them with their name
+NumberOfParameters = M_.param_nbr;
+for ii = 1:NumberOfParameters
+  paramname = M_.param_names{ii};
+  eval([ paramname ' = M_.params(' int2str(ii) ');']);
+end
+
+Omega=(1-alppha)/(1-alppha+alppha*epsilon);
+lambda=(1-theta)*(1-betta*theta)/theta*Omega;      %defined on page 61
+kappa=lambda*(siggma+(varphi+alppha)/(1-alppha));     %defined on page 63
+
+prior_mean=0.03;
+prior_std=0.02;
+log_prior_val=log(normpdf(kappa,prior_mean,prior_std));

tests/kalman/block/fs2000.mod

@@ -0,0 +1,84 @@
+// See fs2000.mod in the examples/ directory for details on the model
+
+var m P c e W R k d n l gy_obs gp_obs y dA;
+varexo e_a e_m;
+
+parameters alp bet gam mst rho psi del;
+
+alp = 0.33;
+bet = 0.99;
+gam = 0.003;
+mst = 1.011;
+rho = 0.7;
+psi = 0.787;
+del = 0.02;
+
+model(block);
+dA = exp(gam+e_a);
+log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
+-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0;
+W = l/n;
+-(psi/(1-psi))*(c*P/(1-n))+l/n = 0;
+R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W;
+1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0;
+c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1);
+P*c = m;
+m-1+d = l;
+e = exp(e_a);
+y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a));
+gy_obs = dA*y/y(-1);
+gp_obs = (P/P(-1))*m(-1)/dA;
+end;
+
+steady_state_model;
+  dA = exp(gam);
+  gst = 1/dA;
+  m = mst;
+  khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1));
+  xist = ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/mst )^(-1);
+  nust = psi*mst^2/( (1-alp)*(1-psi)*bet*gst^alp*khst^alp );
+  n  = xist/(nust+xist);
+  P  = xist + nust;
+  k  = khst*n;
+
+  l  = psi*mst*n/( (1-psi)*(1-n) );
+  c  = mst/P;
+  d  = l - mst + 1;
+  y  = k^alp*n^(1-alp)*gst^alp;
+  R  = mst/bet;
+  W  = l/n;
+  ist  = y-c;
+  q  = 1 - d;
+
+  e = 1;
+  
+  gp_obs = m/dA;
+  gy_obs = dA;
+end;
+
+shocks;
+var e_a; stderr 0.014;
+var e_m; stderr 0.005;
+end;
+
+steady;
+
+check;
+
+estimated_params;
+alp, beta_pdf, 0.356, 0.02;
+bet, beta_pdf, 0.993, 0.002;
+gam, normal_pdf, 0.0085, 0.003;
+mst, normal_pdf, 1.0002, 0.007;
+rho, beta_pdf, 0.129, 0.223;
+psi, beta_pdf, 0.65, 0.05;
+del, beta_pdf, 0.01, 0.005;
+stderr e_a, inv_gamma_pdf, 0.035449, inf;
+stderr e_m, inv_gamma_pdf, 0.008862, inf;
+end;
+
+varobs gp_obs gy_obs;
+
+options_.solve_tolf = 1e-12;
+
+estimation(order=1,datafile='../../fs2000/fsdat_simul',nobs=192,mh_replic=0,mh_nblocks=1,mh_jscale=0.8,consider_only_observed);

tests/kalman/block/fs2000_missing_data.mod

@@ -0,0 +1,83 @@
+// Tests kalman filter with missing observations
+
+var m P c e W R k d n l gy_obs gp_obs y dA;
+varexo e_a e_m;
+
+parameters alp bet gam mst rho psi del;
+
+alp = 0.33;
+bet = 0.99;
+gam = 0.003;
+mst = 1.011;
+rho = 0.7;
+psi = 0.787;
+del = 0.02;
+
+model(block);
+dA = exp(gam+e_a);
+log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
+-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0;
+W = l/n;
+-(psi/(1-psi))*(c*P/(1-n))+l/n = 0;
+R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W;
+1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0;
+c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1);
+P*c = m;
+m-1+d = l;
+e = exp(e_a);
+y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a));
+gy_obs = dA*y/y(-1);
+gp_obs = (P/P(-1))*m(-1)/dA;
+end;
+
+steady_state_model;
+  dA = exp(gam);
+  gst = 1/dA;
+  m = mst;
+  khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1));
+  xist = ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/mst )^(-1);
+  nust = psi*mst^2/( (1-alp)*(1-psi)*bet*gst^alp*khst^alp );
+  n  = xist/(nust+xist);
+  P  = xist + nust;
+  k  = khst*n;
+
+  l  = psi*mst*n/( (1-psi)*(1-n) );
+  c  = mst/P;
+  d  = l - mst + 1;
+  y  = k^alp*n^(1-alp)*gst^alp;
+  R  = mst/bet;
+  W  = l/n;
+  ist  = y-c;
+  q  = 1 - d;
+
+  e = 1;
+  
+  gp_obs = m/dA;
+  gy_obs = dA;
+end;
+
+
+shocks;
+var e_a; stderr 0.014;
+var e_m; stderr 0.005;
+end;
+
+steady;
+
+check;
+
+estimated_params;
+alp, beta_pdf, 0.356, 0.02;
+bet, beta_pdf, 0.993, 0.002;
+gam, normal_pdf, 0.0085, 0.003;
+mst, normal_pdf, 1.0002, 0.007;
+rho, beta_pdf, 0.129, 0.223;
+psi, beta_pdf, 0.65, 0.05;
+del, beta_pdf, 0.01, 0.005;
+stderr e_a, inv_gamma_pdf, 0.035449, inf;
+stderr e_m, inv_gamma_pdf, 0.008862, inf;
+end;
+
+varobs gp_obs gy_obs;
+
+estimation(order=1, datafile='../../fs2000/fsdat_simul_missing_obs', nobs=192, mh_replic=0, mh_nblocks=1, mh_jscale=0.8);

tests/occbin/model_irrcap_twoconstraints/dynrbc_common.inc

@@ -58,9 +58,8 @@ end
 
 % set inputs      
 
-occbin_setup;
-options_.occbin.smoother.debug=1;
-occbin_solver(simul_periods=200,simul_maxit=200,simul_curb_retrench,simul_check_ahead_periods=200);
+occbin_setup(smoother_debug);
+occbin_solver(simul_debug,simul_periodic_solution,simul_periods=200,simul_maxit=200,simul_curb_retrench,simul_check_ahead_periods=200);
 occbin_write_regimes(filename='test',periods=[1:100]);
 
 %% Modify to plot IRFs 

tests/optimal_policy/Ramsey/ramsey_histval.mod

@@ -23,7 +23,7 @@ r=1;
 end;
 
 histval;
-r(0)=1;
+a(0)=-1;
 end;
 
 steady_state_model;
@@ -42,4 +42,5 @@ end;
 options_.dr_display_tol=0;
 planner_objective(ln(c)-phi*((n^(1+gamma))/(1+gamma)));
 ramsey_model(planner_discount=0.99,instruments=(r));
-stoch_simul(order=1,periods=500);
+stoch_simul(order=1);
+evaluate_planner_objective;