Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision
Loading items

Target

Select target project
  • giovanma/dynare
  • giorgiomas/dynare
  • Vermandel/dynare
  • Dynare/dynare
  • normann/dynare
  • MichelJuillard/dynare
  • wmutschl/dynare
  • FerhatMihoubi/dynare
  • sebastien/dynare
  • lnsongxf/dynare
  • rattoma/dynare
  • CIMERS/dynare
  • FredericKarame/dynare
  • SumuduK/dynare
  • MinjeJeon/dynare
  • camilomrch/dynare
  • DoraK/dynare
  • avtishin/dynare
  • selma/dynare
  • claudio_olguin/dynare
  • jeffjiang07/dynare
  • EthanSystem/dynare
  • stepan-a/dynare
  • wjgatt/dynare
  • JohannesPfeifer/dynare
  • gboehl/dynare
  • ebenetce/dynare
  • chskcau/dynare-doc-fixes
28 results
Select Git revision
Loading items
Show changes
Show whitespace changes
Inline Side-by-side
Showing
with 0 additions and 5649 deletions
dynare++/doc/tl.tex → doc/tl.tex
View file @ 654868ce
File moved
dynare++/Makefile.am deleted 100644 → 0
View file @ 53ae5493
SUBDIRS = utils/cc sylv parser/cc tl doc integ kord src tests
dynare++/doc/Makefile.am deleted 100644 → 0
View file @ 53ae5493
EXTRA_DIST = \
dynare++-ramsey.tex \
dynare++-tutorial.tex \
sylvester.tex \
tl.tex \
changelog-old.html \
changelog-sylv-old.html
if HAVE_PDFLATEX
pdf-local: dynare++-ramsey.pdf dynare++-tutorial.pdf sylvester.pdf tl.pdf
endif
%.pdf: %.tex
$(PDFLATEX) $<
$(PDFLATEX) $<
$(PDFLATEX) $<
CLEANFILES = *.pdf *.log *.aux *.out *.toc
dynare++/doc/changelog-old.html deleted 100644 → 0
View file @ 53ae5493
<HTML>
<TITLE>
Dynare++ Change Log
</TITLE>
<!-- $Header$ -->
<BODY>
<TABLE CELLSPACING=2 ALIGN="CENTER" BORDER=1>
<TR>
<TD BGCOLOR="#d0d0d0" WIDTH="85"> <b>Revision</b> </TD>
<TD BGCOLOR="#d0d0d0" WIDTH="85"> <b>Version</b></TD>
<TD BGCOLOR="#d0d0d0" WIDTH="80"> <b>Date</b> </TD>
<TD BGCOLOR="#d0d0d0" WIDTH="600"> <b>Description of changes</b></TD>
</TR>
<TR>
<TD>
<TD>1.3.7
<TD>2008/01/15
<TD>
<TR><TD><TD><TD> <TD> Corrected a serious bug in centralizing a
decision rule. This bug implies that all results based on simulations
of the decision rule were wrong. However results based on stochastic
fix points were correct. Thanks to Wouter J. den Haan and Joris de Wind!
<TR><TD><TD><TD> <TD> Added options --centralize and --no-centralize.
<TR><TD><TD><TD> <TD> Corrected an error of a wrong
variance-covariance matrix in real-time simulations (thanks to Pawel
Zabzcyk).
<TR><TD><TD><TD> <TD> Corrected a bug of integer overflow in refined
faa Di Bruno formula if one of refinements is empty. This bug appeared
when solving models without forward looking variables.
<TR><TD><TD><TD> <TD> Corrected a bug in the Sylvester equation
formerly working only for models with forward looking variables.
<TR><TD><TD><TD> <TD> Corrected a bug in global check printout.
<TR><TD><TD><TD> <TD> Added generating a dump file.
<TR><TD><TD><TD> <TD> Fixed a bug of forgetting repeated assignments
(for example in parameter settings and initval).
<TR><TD><TD><TD> <TD> Added a diff operator to the parser.
<TR>
<TD>1539
<TD>1.3.6
<TD>2008/01/03
<TD>
<TR><TD><TD><TD> <TD> Corrected a bug of segmentation faults for long
names and path names.
<TR><TD><TD><TD> <TD> Changed a way how random numbers are
generated. Dynare++ uses a separate instance of Mersenne twister for
each simulation, this corrects a flaw of additional randomness caused
by operating system scheduler. This also corrects a strange behaviour
of random generator on Windows, where each simulation was getting the
same sequence of random numbers.
<TR><TD><TD><TD> <TD> Added calculation of conditional distributions
controlled by --condper and --condsim.
<TR><TD><TD><TD> <TD> Dropped creating unfoled version of decision
rule at the end. This might consume a lot of memory. However,
simulations might be slower for some models.
<TR>
<TD>1368
<TD>1.3.5
<TD>2007/07/11
<TD>
<TR><TD><TD><TD> <TD> Corrected a bug of useless storing all derivative
indices in a parser. This consumed a lot of memory for large models.
<TR><TD><TD><TD> <TD> Added an option <tt>--ss-tol</tt> controlling a
tolerance used for convergence of a non-linear solver.
<TR><TD><TD><TD> <TD> Corrected buggy interaction of optimal policy
and forward looking variables with more than one period.
<TR><TD><TD><TD> <TD> Variance matrices can be positive
semidefinite. This corrects a bug of throwing an error if estimating
approximation errors on ellipse of the state space with a
deterministic variable.
<TR><TD><TD><TD> <TD> Implemented simulations with statistics
calculated in real-time. Options <tt>--rtsim</tt> and <tt>--rtper</tt>.
<TR>
<TD>1282
<TD>1.3.4
<TD>2007/05/15
<TD>
<TR><TD><TD><TD> <TD>Corrected a bug of wrong representation of NaN in generated M-files.
<TR><TD><TD><TD> <TD>Corrected a bug of occassionaly wrong evaluation of higher order derivatives of integer powers.
<TR><TD><TD><TD> <TD>Implemented automatic handling of terms involving multiple leads.
<TR><TD><TD><TD> <TD>Corrected a bug in the numerical integration, i.e. checking of the precision of the solution.
<TR>
<TD>1090
<TD>1.3.3
<TD>2006/11/20
<TD>
<TR><TD><TD><TD> <TD>Corrected a bug of non-registering an auxiliary variable in initval assignments.
<TR>
<TD>988
<TD>1.3.2
<TD>2006/10/11
<TD>
<TR><TD><TD><TD> <TD>Corrected a few not-serious bugs: segfault on
some exception, error in parsing large files, error in parsing
matrices with comments, a bug in dynare_simul.m
<TR><TD><TD><TD> <TD>Added posibility to specify a list of shocks for
which IRFs are calculated
<TR><TD><TD><TD> <TD>Added --order command line switch
<TR><TD><TD><TD> <TD>Added writing two Matlab files for steady state
calcs
<TR><TD><TD><TD> <TD>Implemented optimal policy using keyword
planner_objective and planner_discount
<TR><TD><TD><TD> <TD>Implemented an R interface to Dynare++ algorithms
(Tamas Papp)
<TR><TD><TD><TD> <TD>Highlevel code reengineered to allow for
different model inputs
<TR>
<TD>799
<TD>1.3.1
<TD>2006/06/13
<TD>
<TR><TD><TD><TD> <TD>Corrected few bugs: in error functions, in linear algebra module.
<TR><TD><TD><TD> <TD>Updated dynare_simul.
<TR><TD><TD><TD> <TD>Updated the tutorial.
<TR><TD><TD><TD> <TD>Corrected an error in summing up tensors where
setting up the decision rule derivatives. Thanks to Michel
Juillard. The previous version was making deterministic effects of
future volatility smaller than they should be.
<TR>
<TD>766
<TD>1.3.0
<TD>2006/05/22
<TD>
<TR><TD><TD><TD> <TD>The non-linear solver replaced with a new one.
<TR><TD><TD><TD> <TD>The parser and derivator replaced with a new
code. Now it is possible to put expressions in parameters and initval
sections.
<TR>
<TD>752
<TD>1.2.2
<TD>2006/05/22
<TD>
<TR><TD><TD><TD> <TD>Added an option triggering/suppressing IRF calcs..
<TR><TD><TD><TD> <TD>Newton algortihm is now used for fix-point calculations.
<TR><TD><TD><TD> <TD> Vertical narrowing of tensors in Faa Di Bruno
formula to avoid multiplication with zeros..
<TR>
<TD>436
<TD>1.2.1
<TD>2005/08/17
<TD>
<TR><TD><TD><TD> <TD>Faa Di Bruno for sparse matrices optimized. The
implementation now accommodates vertical refinement of function stack
in order to fit a corresponding slice to available memory. In
addition, zero slices are identified. For some problems, this implies
significant speedup.
<TR><TD><TD><TD> <TD>Analytic derivator speedup.
<TR><TD><TD><TD> <TD>Corrected a bug in the threading code. The bug
stayed concealed in Linux 2.4.* kernels, and exhibited in Linux 2.6.*,
which has a different scheduling. This correction also allows using
detached threads on Windows.
<TR>
<TD>410
<TD>1.2
<TD>2005/07/29
<TD>
<TR><TD><TD><TD> <TD>Added Dynare++ tutorial.
<TR><TD><TD><TD> <TD>Changed and enriched contents of MAT-4 output
file.
<TR><TD><TD><TD> <TD>Corrected a bug of wrong variable indexation
resulting in an exception. The error occurred if a variable appeared
at time t-1 or t+1 and not at t.
<TR><TD><TD><TD> <TD>Added Matlab interface, which allows simulation
of a decision rule in Matlab.
<TR><TD><TD><TD> <TD>Got rid of Matrix Template Library.
<TR><TD><TD><TD> <TD>Added checking of model residuals by the
numerical integration. Three methods: checking along simulation path,
checking along shocks, and on ellipse of states.
<TR><TD><TD><TD> <TD>Corrected a bug in calculation of higher moments
of Normal dist.
<TR><TD><TD><TD> <TD>Corrected a bug of wrong drawing from Normal dist
with non-zero covariances.
<TR><TD><TD><TD>
<TD>Added numerical integration module. Product and Smolyak
quadratures over Gauss-Hermite and Gauss-Legendre, and quasi Monte
Carlo.
<TR>
<TD>152
<TD>1.1
<TD>2005/04/22
<TD>
<TR><TD><TD><TD>
<TD>Added a calculation of approximation at a stochastic steady state
(still experimental).
<TR><TD><TD><TD>
<TD>Corrected a bug in Cholesky decomposition of variance-covariance
matrix with off-diagonal elements.
<TR>
<TD>89
<TD>1.01
<TD>2005/02/23
<TD>
<TR><TD><TD><TD>
<TD>Added version printout.
<TR><TD><TD><TD>
<TD>Corrected the bug of multithreading support for P4 HT processors running on Win32.
<TR><TD><TD><TD>
<TD>Enhanced Kronecker product code resulting in approx. 20% speedup.
<TR><TD><TD><TD>
<TD>Implemented vertical stack container refinement, and another
method for sparse folded Faa Di Bruno (both not used yet).
<TR>
<TD>5
<TD>1.0
<TD>2005/02/23
<TD>The first released version.
</TABLE>
</BODY>
</HTML>
dynare++/doc/changelog-sylv-old.html deleted 100644 → 0
View file @ 53ae5493
<HTML>
<!-- $Header: /var/lib/cvs/dynare_cpp/sylv/change_log.html,v 1.1.1.1 2004/06/04 13:00:05 kamenik Exp $ -->
<!-- Tag $Name: $ -->
<TITLE>
Sylvester Solver Change Log
</TITLE>
<BODY>
<TABLE CELLSPACING=2 ALIGN="CENTER" BORDER=1>
<TR>
<TD BGCOLOR="#d0d0d0" WIDTH="85"> Tag </TD>
<TD BGCOLOR="#d0d0d0" WIDTH="80"> Date </TD>
<TD BGCOLOR="#d0d0d0" WIDTH="600"> Description/Changes</TD>
</TR>
<TR>
<TD></TD>
<TD>2003/09/10</TD>
<TD>Initial version solving triangular system put to repository</TD>
</TR>
<TR>
<TD></TD>
<TD></TD>
<TD>Implemented solution of general case.</TD>
</TR>
<TR>
<TD></TD>
<TD></TD>
<TD>Implemented a memory pool (Paris).</TD>
</TR>
<TR>
<TD></TD>
<TD></TD>
<TD>Implemented MEX interface to the routine (Paris).</TD>
</TR>
<TR>
<TD></TD>
<TD></TD>
<TD>Implemented QuasiTriangularZero (Paris) (not fully used yet).</TD>
</TR>
<TR>
<TD>rel-1</TD>
<TD>2003/10-02</TD>
<TD>Version sent to Michel.</TD>
</TR>
<TR>
<TD></TD>
<TD></TD>
<TD>Inheritance streamlined, QuasiTriangular inherits from GeneralMatrix.</TD>
</TR>
<TR>
<TD></TD>
<TD></TD>
<TD>Implemented block diagonalization algorithm.</TD>
</TR>
<TR>
<TD></TD>
<TD></TD>
<TD>Solution routines rewritten so that the output rewrites input,
considerable memory improvement.</TD>
</TR>
<TR>
<TD></TD>
<TD></TD>
<TD>MEX interface now links with LAPACK library from Matlab.</TD>
</TR>
<TR>
<TD></TD>
<TD></TD>
<TD>Added a hack to MEX library loading in order to avoid Matlab crash in Wins.</TD>
</TR>
<TR>
<TD>rel-2</TD>
<TD>2003/10/15</TD>
<TD>Version sent to Michel.</TD>
</TR>
<TR>
<TD></TD>
<TD></TD>
<TD>KronUtils now rewrite input by output using less memory.</TD>
</TR>
<TR>
<TD></TD>
<TD></TD>
<TD>Added iterative solution algorithm (doubling).</TD>
</TR>
<TR>
<TD></TD>
<TD></TD>
<TD>Introduced abstraction for set of parameters (SylvParams).</TD>
</TR>
<TR>
<TD></TD>
<TD></TD>
<TD>Algorithm enabled to solve problems with singular C.</TD>
</TR>
<TR>
<TD></TD>
<TD></TD>
<TD>Implemented a class chooser chossing between QuasiTriangularZero,
and QuasiTriangular (padded with zero) depending on size of the
problem. Full use of QuasiTriangularZero.</TD>
</TR>
<TR>
<TD></TD>
<TD></TD>
<TD>Reimplemented QuasiTriangular::solve, offdiagonal elements are
eleiminated by gauss with partial pivoting, not by transformation of
complex eigenvalues. More stable for ill conditioned eigenvalues.</TD>
</TR>
<TR>
<TD></TD>
<TD></TD>
<TD>Reimplemented calculation of eliminating vectors, much more
numerically stable now.</TD>
</TR>
<TR>
<TD></TD>
<TD></TD>
<TD>Implemented algorithm for ordering of eigenvalues (not used now,
no numerical improvements).</TD>
</TR>
<TR>
<TD>rel-3</TD>
<TD>2003/12/4</TD>
<TD>Version sent to Michel.</TD>
</TR>
<TR>
<TD></TD>
<TD></TD>
<TD>GeneralMatrix separated for use outside, in sylv module we use
its subclass SylvMatrix. Implemented ConstGeneralMatrix (useful outside).
</TD>
</TR>
<TR>
<TD>rel-4</TD>
<TD>2004/6/4</TD>
<TD>Version, which was moved to pythie.cepremap.cnrs.fr repository.</TD>
</TR>
</TABLE>
</BODY>
</HTML>
dynare++/doc/dynare++-ramsey.tex deleted 100644 → 0
View file @ 53ae5493
\documentclass[10pt]{article}
\usepackage{array,natbib,times}
\usepackage{amsmath, amsthm, amssymb}
%\usepackage[pdftex,colorlinks]{hyperref}
\begin{document}
\title{Implementation of Ramsey Optimal Policy in Dynare++, Timeless Perspective}
\author{Ondra Kamen\'\i k}
\date{June 2006}
\maketitle
\textbf{Abstract:} This document provides a derivation of Ramsey
optimal policy from timeless perspective and describes its
implementation in Dynare++.
\section{Derivation of the First Order Conditions}
Let us start with an economy populated by agents who take a number of
variables exogenously, or given. These may include taxes or interest
rates for example. These variables can be understood as decision (or control)
variables of the timeless Ramsey policy (or social planner). The agent's
information set at time $t$ includes mass-point distributions of these
variables for all times after $t$. If $i_t$ denotes an interest rate
for example, then the information set $I_t$ includes
$i_{t|t},i_{t+1|t},\ldots,i_{t+k|t},\ldots$ as numbers. In addition
the information set includes all realizations of past exogenous
innovations $u_\tau$ for $\tau=t,t-1,\ldots$ and distibutions
$u_\tau\sim N(0,\Sigma)$ for $\tau=t+1,\ldots$. These information sets will be denoted $I_t$.
An information set including only the information on past realizations
of $u_\tau$ and future distributions of $u_\tau\sim N(0\sigma)$ will
be denoted $J_t$. We will use the following notation for expectations
through these sets:
\begin{eqnarray*}
E^I_t[X] &=& E(X|I_t)\\
E^J_t[X] &=& E(X|J_t)
\end{eqnarray*}
The agents optimize taking the decision variables of the social
planner at $t$ and future as given. This means that all expectations
they form are conditioned on the set $I_t$. Let $y_t$ denote a vector
of all endogenous variables including the planer's decision
variables. Let the number of endogenous variables be $n$. The economy
can be described by $m$ equations including the first order conditions
and transition equations:
\begin{equation}\label{constr}
E_t^I\left[f(y_{t-1},y_t,y_{t+1},u_t)\right] = 0.
\end{equation}
This lefts $n-m$
the planner's control variables. The solution of this problem is a
decision rule of the form:
\begin{equation}\label{agent_dr}
y_t=g(y_{t-1},u_t,c_{t|t},c_{t+1|t},\ldots,c_{t+k|t},\ldots),
\end{equation}
where $c$ is a vector of planner's control variables.
Each period the social planner chooses the vector $c_t$ to maximize
his objective such that \eqref{agent_dr} holds for all times following
$t$. This would lead to $n-m$ first order conditions with respect to
$c_t$. These first order conditions would contain unknown derivatives
of endogenous variables with respect to $c$, which would have to be
retrieved from the implicit constraints \eqref{constr} since the
explicit form \eqref{agent_dr} is not known.
The other way to proceed is to assume that the planner is so dumb that
he is not sure what are his control variables. So he optimizes with
respect to all $y_t$ given the constraints \eqref{constr}. If the
planner's objective is $b(y_{t-1},y_t,y_{t+1},u_t)$ with a discount rate
$\beta$, then the optimization problem looks as follows:
\begin{align}
\max_{\left\{y_\tau\right\}^\infty_t}&E_t^J
\left[\sum_{\tau=t}^\infty\beta^{\tau-t}b(y_{\tau-1},y_\tau,y_{\tau+1},u_\tau)\right]\notag\\
&\rm{s.t.}\label{planner_optim}\\
&\hskip1cm E^I_\tau\left[f(y_{\tau-1},y_\tau,y_{\tau+1},u_\tau)\right]=0\quad\rm{for\ }
\tau=\ldots,t-1,t,t+1,\ldots\notag
\end{align}
Note two things: First, each constraint \eqref{constr} in
\eqref{planner_optim} is conditioned on $I_\tau$ not $I_t$. This is
very important, since the behaviour of agents at period $\tau=t+k$ is
governed by the constraint using expectations conditioned on $t+k$,
not $t$. The social planner knows that at $t+k$ the agents will use
all information available at $t+k$. Second, the constraints for the
planner's decision made at $t$ include also constraints for agent's
behaviour prior to $t$. This is because the agent's decision rules are
given in the implicit form \eqref{constr} and not in the explicit form
\eqref{agent_dr}.
Using Lagrange multipliers, this can be rewritten as
\begin{align}
\max_{y_t}E_t^J&\left[\sum_{\tau=t}^\infty\beta^{\tau-t}b(y_{\tau-1},y_\tau,y_{\tau+1},u_\tau)\right.\notag\\
&\left.+\sum_{\tau=-\infty}^{\infty}\beta^{\tau-t}\lambda^T_\tau E_\tau^I\left[f(y_{\tau-1},y_\tau,y_{\tau+1},u_\tau)\right]\right],
\label{planner_optim_l}
\end{align}
where $\lambda_t$ is a vector of Lagrange multipliers corresponding to
constraints \eqref{constr}. Note that the multipliers are multiplied
by powers of $\beta$ in order to make them stationary. Taking a
derivative wrt $y_t$ and putting it to zero yields the first order
conditions of the planner's problem:
\begin{align}
E^J_t\left[\vphantom{\frac{\int^(_)}{\int^(\_)}}\right.&\frac{\partial}{\partial y_t}b(y_{t-1},y_t,y_{t+1},u_t)+
\beta L^{+1}\frac{\partial}{\partial y_{t-1}}b(y_{t-1},y_t,y_{t+1},u_t)\notag\\
&+\beta^{-1}\lambda_{t-1}^TE^I_{t-1}\left[L^{-1}\frac{\partial}{\partial y_{t+1}}f(y_{t-1},y_t,y_{t+1},u_t)\right]\notag\\
&+\lambda_t^TE^I_t\left[\frac{\partial}{\partial y_{t}}f(y_{t-1},y_t,y_{t+1},u_t)\right]\notag\\
&+\beta\lambda_{t+1}^TE^I_{t+1}\left[L^{+1}\frac{\partial}{\partial y_{t-1}}f(y_{t-1},y_t,y_{t+1},u_t)\right]
\left.\vphantom{\frac{\int^(_)}{\int^(\_)}}\right]
= 0,\label{planner_optim_foc}
\end{align}
where $L^{+1}$ and $L^{-1}$ are one period lead and lag operators respectively.
Now we have to make a few assertions concerning expectations
conditioned on the different information sets to simplify
\eqref{planner_optim_foc}. Recall the formula for integration through
information on which another expectation is conditioned, this is:
$$E\left[E\left[u|v\right]\right] = E[u],$$
where the outer expectation integrates through $v$. Since $J_t\subset
I_t$, by easy application of the above formula we obtain
\begin{eqnarray}
E^J_t\left[E^I_t\left[X\right]\right] &=& E^J_t\left[X\right]\quad\rm{and}\notag\\
E^J_t\left[E^I_{t-1}\left[X\right]\right] &=& E^J_t\left[X\right]\label{e_iden}\\
E^J_t\left[E^I_{t+1}\left[X\right]\right] &=& E^J_{t+1}\left[X\right]\notag
\end{eqnarray}
Now, the last term of \eqref{planner_optim_foc} needs a special
attention. It is equal to
$E^J_t\left[\beta\lambda^T_{t+1}E^I_{t+1}[X]\right]$. If we assume
that the problem \eqref{planner_optim} has a solution, then there is a
deterministic function from $J_{t+1}$ to $\lambda_{t+1}$ and so
$\lambda_{t+1}\in J_{t+1}\subset I_{t+1}$. And the last term is equal
to $E^J_{t}\left[E^I_{t+1}[\beta\lambda^T_{t+1}X]\right]$, which is
$E^J_{t+1}\left[\beta\lambda^T_{t+1}X\right]$. This term can be
equivalently written as
$E^J_{t}\left[\beta\lambda^T_{t+1}E^J_{t+1}[X]\right]$. The reason why
we write the term in this way will be clear later. All in all, we have
\begin{align}
E^J_t\left[\vphantom{\frac{\int^(_)}{\int^(\_)}}\right.&\frac{\partial}{\partial y_t}b(y_{t-1},y_t,y_{t+1},u_t)+
\beta L^{+1}\frac{\partial}{\partial y_{t-1}}b(y_{t-1},y_t,y_{t+1},u_t)\notag\\
&+\beta^{-1}\lambda_{t-1}^TL^{-1}\frac{\partial}{\partial y_{t+1}}f(y_{t-1},y_t,y_{t+1},u_t)\notag\\
&+\lambda_t^T\frac{\partial}{\partial y_{t}}f(y_{t-1},y_t,y_{t+1},u_t)\notag\\
&+\beta\lambda_{t+1}^TE^J_{t+1}\left[L^{+1}\frac{\partial}{\partial y_{t-1}}f(y_{t-1},y_t,y_{t+1},u_t)\right]
\left.\vphantom{\frac{\int^(_)}{\int^(\_)}}\right]
= 0.\label{planner_optim_foc2}
\end{align}
Note that we have not proved that \eqref{planner_optim_foc} and
\eqref{planner_optim_foc2} are equivalent. We proved only that if
\eqref{planner_optim_foc} has a solution, then
\eqref{planner_optim_foc2} is equivalent (and has the same solution).
%%- \section{Implementation}
%%-
%%- The user inputs $b(y_{t-1},y_t,y_{t+1},u_t)$, $\beta$, and agent's
%%- first order conditions \eqref{constr}. The algorithm has to produce
%%- \eqref{planner_optim_foc2}.
%%-
\end{document}
dynare++/doc/dynare++-tutorial.tex deleted 100644 → 0
View file @ 53ae5493
\documentclass[10pt]{article}
\usepackage{array,natbib}
\usepackage{amsmath, amsthm, amssymb}
\usepackage[pdftex,colorlinks]{hyperref}
\begin{document}
\title{DSGE Models with Dynare++. A Tutorial.}
\author{Ondra Kamen\'\i k}
\date{February 2011, updated August 2016}
\maketitle
\tableofcontents
\section{Setup}
The Dynare++ setup procedure is pretty straightforward as Dynare++ is included in the Dynare installation
packages which can be downloaded from \url{http://www.dynare.org}. Take the following steps:
\begin{enumerate}
\item Add the {\tt dynare++} subdirectory of the root Dynare installation directory to the your
operating system path. This ensures that your OS will find the {\tt dynare++} executable.
\item If you have MATLAB and want to run custom simulations (see \ref{custom}),
then you need to add to your MATLAB path the {\tt dynare++} subdirectory of
the root Dynare installation directory, and also directory containing the
\texttt{dynare\_simul\_} MEX file (note the trailing underscore). The easiest
way to add the latter is to run Dynare once in your MATLAB session (even
without giving it any MOD file).
\end{enumerate}
\section{Sample Session}
As an example, let us take a simple DSGE model with time to build, whose dynamic
equilibrium is described by the following first order conditions:
\begin{align*}
&c_t\theta h_t^{1+\psi} = (1-\alpha)y_t\cr
&\beta E_t\left[\frac{\exp(b_t)c_t}{\exp(b_{t+1})c_{t+1}}
\left(\exp(b_{t+1})\alpha\frac{y_{t+1}}{k_{t+1}}+1-\delta\right)\right]=1\cr
&y_t=\exp(a_t)k_t^\alpha h_t^{1-\alpha}\cr
&k_{t}=\exp(b_{t-1})(y_{t-1}-c_{t-1})+(1-\delta)k_{t-1}\cr
&a_t=\rho a_{t-1}+\tau b_{t-1}+\epsilon_t\cr
&b_t=\tau a_{t-1}+\rho b_{t-1}+\nu_t
\end{align*}
\label{timing}
The convention is that the timing of a variable reflects when this variable
is decided. Dynare++ therefore uses a ``stock at the end of the
period'' notation for predetermined state variables (see the Dynare manual for details).
The timing of this model is that the exogenous shocks $\epsilon_t$,
and $\nu_t$ are observed by agents in the beginning of period $t$ and
before the end of period $t$ all endogenous variables with index $t$
are decided. The expectation operator $E_t$ works over the information
accumulated just before the end of the period $t$ (this includes
$\epsilon_t$, $\nu_t$ and all endogenous variables with index $t$).
The exogenous shocks $\epsilon_t$ and $\nu_t$ are supposed to be
serially uncorrelated with zero means and time-invariant
variance-covariance matrix. In Dynare++, these variables are called
exogenous; all other variables are endogenous. Now we are prepared to
start writing a model file for Dynare++, which is an ordinary text
file and could be created with any text editor.
The model file starts with a preamble declaring endogenous and
exogenous variables, parameters, and setting values of the
parameters. Note that one can put expression on right hand sides. The
preamble follows:
{\small
\begin{verbatim}
var Y, C, K, A, H, B;
varexo EPS, NU;
parameters beta, rho, alpha, delta, theta, psi, tau;
alpha = 0.36;
rho = 0.95;
tau = 0.025;
beta = 1/(1.03^0.25);
delta = 0.025;
psi = 0;
theta = 2.95;
\end{verbatim}
}
The section setting values of the parameters is terminated by a
beginning of the {\tt model} section, which states all the dynamic
equations. A timing convention of a Dynare++ model is the same as the
timing of our example model, so we may proceed with writing the model
equations. The time indexes of $c_{t-1}$, $c_t$, and $c_{t+1}$ are
written as {\tt C(-1)}, {\tt C}, and {\tt C(1)} resp. The {\tt model}
section looks as follows:
{\small
\begin{verbatim}
model;
C*theta*H^(1+psi) = (1-alpha)*Y;
beta*exp(B)*C/exp(B(1))/C(1)*
(exp(B(1))*alpha*Y(1)/K(1)+1-delta) = 1;
Y = exp(A)*K^alpha*H^(1-alpha);
K = exp(B(-1))*(Y(-1)-C(-1)) + (1-delta)*K(-1);
A = rho*A(-1) + tau*B(-1) + EPS;
B = tau*A(-1) + rho*B(-1) + NU;
end;
\end{verbatim}
}
At this point, almost all information that Dynare++ needs has been
provided. Only three things remain to be specified: initial values of
endogenous variables for non-linear solver, variance-covariance matrix
of the exogenous shocks and order of the Taylor approximation. Since
the model is very simple, there is a closed form solution for the
deterministic steady state. We use it as initial values for the
non-linear solver. Note that the expressions on the right hand-sides in
{\tt initval} section can reference values previously calculated. The
remaining portion of the model file looks as follows:
{\small
\begin{verbatim}
initval;
A = 0;
B = 0;
H = ((1-alpha)/(theta*(1-(delta*alpha)
/(1/beta-1+delta))))^(1/(1+psi));
Y = (alpha/(1/beta-1+delta))^(alpha/(1-alpha))*H;
K = alpha/(1/beta-1+delta)*Y;
C = Y - delta*K;
end;
vcov = [
0.0002 0.00005;
0.00005 0.0001
];
order = 7;
\end{verbatim}
}
Note that the order of rows/columns of the variance-covariance matrix
corresponds to the ordering of exogenous variables in the {\tt varexo}
declaration. Since the {\tt EPS} was declared first, its variance is
$0.0002$, and the variance of {\tt NU} is $0.0001$.
Let the model file be saved as {\tt example1.mod}. Now we are prepared
to solve the model. At the operating system command
prompt\footnote{Under Windows it is a {\tt cmd} program, under Unix it
is any shell} we issue a command:
{\small
\begin{verbatim}
dynare++ example1.mod
\end{verbatim}
}
When the program is finished, it produces two output files: a journal
file {\tt example1.jnl} and a Matlab MAT-4 {\tt example1.mat}. The
journal file contains information about time, memory and processor
resources needed for all steps of solution. The output file is more
interesting. It contains various simulation results. It can be loaded
into Matlab or Scilab and examined.%
\footnote{For Matlab {\tt load example1.mat}, for Scilab {\tt
mtlb\_load example1.mat}} The following examples are done in Matlab,
everything would be very similar in Scilab.
Let us first examine the contents of the MAT file:
{\small
\begin{verbatim}
>> load example1.mat
>> who
Your variables are:
dyn_g_1 dyn_i_Y dyn_npred
dyn_g_2 dyn_irfm_EPS_mean dyn_nstat
dyn_g_3 dyn_irfm_EPS_var dyn_shocks
dyn_g_4 dyn_irfm_NU_mean dyn_ss
dyn_g_5 dyn_irfm_NU_var dyn_state_vars
dyn_i_A dyn_irfp_EPS_mean dyn_steady_states
dyn_i_B dyn_irfp_EPS_var dyn_vars
dyn_i_C dyn_irfp_NU_mean dyn_vcov
dyn_i_EPS dyn_irfp_NU_var dyn_vcov_exo
dyn_i_H dyn_mean
dyn_i_K dyn_nboth
dyn_i_NU dyn_nforw
\end{verbatim}
}
All the variables coming from one MAT file have a common prefix. In
this case it is {\tt dyn}, which is Dynare++ default. The prefix can
be changed, so that the multiple results could be loaded into one Matlab
session.
In the default setup, Dynare++ solves the Taylor approximation to the
decision rule and calculates unconditional mean and covariance of the
endogenous variables, and generates impulse response functions. The
mean and covariance are stored in {\tt dyn\_mean} and {\tt
dyn\_vcov}. The ordering of the endogenous variables is given by {\tt
dyn\_vars}.
In our example, the ordering is
{\small
\begin{verbatim}
>> dyn_vars
dyn_vars =
H
A
Y
C
K
B
\end{verbatim}
}
and unconditional mean and covariance are
{\small
\begin{verbatim}
>> dyn_mean
dyn_mean =
0.2924
0.0019
1.0930
0.8095
11.2549
0.0011
>> dyn_vcov
dyn_vcov =
0.0003 0.0006 0.0016 0.0004 0.0060 0.0004
0.0006 0.0024 0.0059 0.0026 0.0504 0.0012
0.0016 0.0059 0.0155 0.0069 0.1438 0.0037
0.0004 0.0026 0.0069 0.0040 0.0896 0.0016
0.0060 0.0504 0.1438 0.0896 2.1209 0.0405
0.0004 0.0012 0.0037 0.0016 0.0405 0.0014
\end{verbatim}
}
The ordering of the variables is also given by indexes starting with
{\tt dyn\_i\_}. Thus the mean of capital can be retrieved as
{\small
\begin{verbatim}
>> dyn_mean(dyn_i_K)
ans =
11.2549
\end{verbatim}
}
\noindent and covariance of labor and capital by
{\small
\begin{verbatim}
>> dyn_vcov(dyn_i_K,dyn_i_H)
ans =
0.0060
\end{verbatim}
}
The impulse response functions are stored in matrices as follows
\begin{center}
\begin{tabular}{|l|l|}
\hline
matrix& response to\\
\hline
{\tt dyn\_irfp\_EPS\_mean}& positive impulse to {\tt EPS}\\
{\tt dyn\_irfm\_EPS\_mean}& negative impulse to {\tt EPS}\\
{\tt dyn\_irfp\_NU\_mean}& positive impulse to {\tt NU}\\
{\tt dyn\_irfm\_NU\_mean}& negative impulse to {\tt NU}\\
\hline
\end{tabular}
\end{center}
All shocks sizes are one standard error. Rows of the matrices
correspond to endogenous variables, columns correspond to
periods. Thus capital response to a positive shock to {\tt EPS} can be
plotted as
{\small
\begin{verbatim}
plot(dyn_irfp_EPS_mean(dyn_i_K,:));
\end{verbatim}
}
The data is in units of the respective variables, so in order to plot
the capital response in percentage changes from the decision rule's
fix point (which is a vector {\tt dyn\_ss}), one has to issue the
commands:
{\small
\begin{verbatim}
Kss=dyn_ss(dyn_i_K);
plot(100*dyn_irfp_EPS_mean(dyn_i_K,:)/Kss);
\end{verbatim}
}
The plotted impulse response shows that the model is pretty persistent
and that the Dynare++ default for a number of simulated periods is not
sufficient. In addition, the model persistence puts in doubt also a
number of simulations. The Dynare++ defaults can be changed when
calling Dynare++, in operating system's command prompt, we issue a
command:
{\small
\begin{verbatim}
dynare++ --per 300 --sim 150 example1.mod
\end{verbatim}
}
\noindent This sets the number of simulations to $150$ and the number
of periods to $300$ for each simulation giving $45000$ total simulated
periods.
\section{Sample Optimal Policy Session}
\label{optim_tut}
Suppose that one wants to solve the following optimal policy problem
with timeless perspective.\footnote{See \ref{ramsey} on how to solve
Ramsey optimality problem within this framework} The following
optimization problem is how to choose capital taxes financing public
good to maximize agent's utility from consumption good and public
good. The problem takes the form:
\begin{align*}
\max_{\{\tau_t\}_{t_0}^\infty}
E_{t_0}\sum_{t=t_0}^\infty &\beta^{t-t_0}\left(u(c_t)+av(g_t)\right)\\
\hbox{subject\ to}&\\
u'(c_t) &=
\beta E_t\left[u'(c_{t+1})\left(1-\delta+f'(k_{t+1})(1-\alpha\tau_{t+1})\right)\right]\\
K_t &= (1-\delta)K_{t-1} + (f(K_{t-1}) - c_{t-1} - g_{t-1})\\
g_t &= \tau_t\alpha f(K_t),\\
\hbox{where\ } t & = \ldots,t_0-1,t_0,t_0+1,\ldots
\end{align*}
$u(c_t)$ is utility from consuming the consumption good, $v(g_t)$ is
utility from consuming the public good, $f(K_t)$ is a production
function $f(K_t) = Z_tK_t^\alpha$. $Z_t$ is a technology shock modeled
as AR(1) process. The three constraints come from the first order
conditions of a representative agent. We suppose that it pursues a
different objective, namely lifetime utility involving only
consumption $c_t$. The representative agents chooses between
consumption and investment. It rents the capital to firms and supplies
constant amount of labour. All output is paid back to consumer in form
of wage and capital rent. Only the latter is taxed. We suppose that
the optimal choice has been taking place from infinite past and will
be taking place for ever. Further we suppose the same about the
constraints.
Let us choose the following functional forms:
\begin{eqnarray*}
u(c_t) &=& \frac{c_t^{1-\eta}}{1-\eta}\\
v(g_t) &=& \frac{g_t^{1-\phi}}{1-\phi}\\
f(K_t) &=& K_t^\alpha
\end{eqnarray*}
Then the problem can be coded into Dynare++ as follows. We start with
a preamble which states all the variables, shocks and parameters:
{\small
\begin{verbatim}
var C G K TAU Z;
varexo EPS;
parameters eta beta alpha delta phi a rho;
eta = 2;
beta = 0.99;
alpha = 0.3;
delta = 0.10;
phi = 2.5;
a = 0.1;
rho = 0.7;
\end{verbatim}
}
Then we specify the planner's objective and the discount factor in the
objective. The objective is an expression (possibly including also
variable leads and lags), and the discount factor must be one single
declared parameter:
{\small
\begin{verbatim}
planner_objective C^(1-eta)/(1-eta) + a*G^(1-phi)/(1-phi);
planner_discount beta;
\end{verbatim}
}
The model section will contain only the constraints of the social
planner. These are capital accumulation, identity for the public
product, AR(1) process for $Z_t$ and the first order condition of the
representative agent (with different objective).
{\small
\begin{verbatim}
model;
K = (1-delta)*K(-1) + (exp(Z(-1))*K(-1)^alpha - C(-1) - G(-1));
G = TAU*alpha*K^alpha;
Z = rho*Z(-1) + EPS;
C^(-eta) = beta*C(+1)^(-eta)*(1-delta +
exp(Z(+1))*alpha*K(+1)^(alpha-1)*(1-alpha*TAU(+1)));
end;
\end{verbatim}
}
Now we have to provide a good guess for non-linear solver calculating
the deterministic steady state. The model's steady state has a closed
form solution if the taxes are known. So we provide a guess for
taxation {\tt TAU} and then use the closed form solution for capital,
public good and consumption:\footnote{Initial guess for Lagrange
multipliers and some auxiliary variables is calculated automatically. See
\ref{opt_init} for more details.}
{\small
\begin{verbatim}
initval;
TAU = 0.70;
K = ((delta+1/beta-1)/(alpha*(1-alpha*TAU)))^(1/(alpha-1));
G = TAU*alpha*K^alpha;
C = K^alpha - delta*K - G;
Z = 0;
\end{verbatim}
}
Finally, we have to provide the order of approximation, and the
variance-covariance matrix of the shocks (in our case we have only one
shock):
{\small
\begin{verbatim}
order = 4;
vcov = [
0.01
];
\end{verbatim}
}
After this model file has been run, we can load the resulting MAT-file
into the Matlab (or Scilab) and examine its contents:
{\small
\begin{verbatim}
>> load kp1980_2.mat
>> who
Your variables are:
dyn_g_1 dyn_i_MULT1 dyn_nforw
dyn_g_2 dyn_i_MULT2 dyn_npred
dyn_g_3 dyn_i_MULT3 dyn_nstat
dyn_g_4 dyn_i_TAU dyn_shocks
dyn_i_AUX_3_0_1 dyn_i_Z dyn_ss
dyn_i_AUX_4_0_1 dyn_irfm_EPS_mean dyn_state_vars
dyn_i_C dyn_irfm_EPS_var dyn_steady_states
dyn_i_EPS dyn_irfp_EPS_mean dyn_vars
dyn_i_G dyn_irfp_EPS_var dyn_vcov
dyn_i_K dyn_mean dyn_vcov_exo
dyn_i_MULT0 dyn_nboth
\end{verbatim}
}
The data dumped into the MAT-file have the same structure as in the
previous example of this tutorial. The only difference is that
Dynare++ added a few more variables. Indeed:
{\small
\begin{verbatim}
>> dyn_vars
dyn_vars =
MULT1
G
MULT3
C
K
Z
TAU
AUX_3_0_1
AUX_4_0_1
MULT0
MULT2
\end{verbatim}
}
Besides the five variables declared in the model ({\tt C}, {\tt G},
{\tt K}, {\tt TAU}, and {\tt Z}), Dy\-na\-re++ added 6 more, four as Lagrange
multipliers of the four constraints, two as auxiliary variables for
shifting in time. See \ref{aux_var} for more details.
The structure and the logic of the MAT-file is the same as these new 6
variables were declared in the model file and the file is examined in
the same way.
For instance, let us examine the Lagrange multiplier of the optimal
policy associated with the consumption first order condition. Recall
that the consumers' objective is different from the policy
objective. Therefore, the constraint will be binding and the
multiplier will be non-zero. Indeed, its deterministic steady state,
fix point and mean are as follows:
{\small
\begin{verbatim}
>> dyn_steady_states(dyn_i_MULT3,1)
ans =
-1.3400
>> dyn_ss(dyn_i_MULT3)
ans =
-1.3035
>> dyn_mean(dyn_i_MULT3)
ans =
-1.3422
\end{verbatim}
}
\section{What Dynare++ Calculates}
\label{dynpp_calc}
Dynare++ solves first order conditions of a DSGE model in the recursive form:
\begin{equation}\label{focs}
E_t[f(y^{**}_{t+1},y_t,y^*_{t-1},u_t)]=0,
\end{equation}
where $y$ is a vector of endogenous variables, and $u$ a vector of
exogenous variables. Some of elements of $y$ can occur at time $t+1$,
these are $y^{**}$. Elements of $y$ occurring at time $t-1$ are denoted
$y^*$. The exogenous shocks are supposed to be serially independent
and normally distributed $u_t\sim N(0,\Sigma)$.
The solution of this dynamic system is a decision rule
\[
y_t=g(y^*_{t-1},u_t)
\]
Dynare++ calculates a Taylor approximation of this decision rule of a
given order. The approximation takes into account deterministic
effects of future volatility, so a point about which the Taylor
approximation is done will be different from the fix point $y$ of the rule
yielding $y=g(y^*,0)$.
The fix point of a rule corresponding to a model with $\Sigma=0$ is
called {\it deterministic steady state} denoted as $\bar y$. In
contrast to deterministic steady state, there is no consensus in
literature how to call a fix point of the rule corresponding to a
model with non-zero $\Sigma$. I am tempted to call it {\it stochastic
steady state}, however, it might be confused with unconditional mean
or with steady distribution. So I will use a term {\it fix point} to
avoid a confusion.
By default, Dynare++ solves the Taylor approximation about the
deterministic steady state. Alternatively, Dynare++ can split the
uncertainty to a few steps and take smaller steps when calculating the
fix points. This is controlled by an option {\tt --steps}. For the
brief description of the second method, see \ref{multistep_alg}.
\subsection{Decision Rule Form}
\label{dr_form}
In case of default solution algorithm (approximation about the
deterministic steady state $\bar y$), Dynare++ calculates the higher
order derivatives of the equilibrium rule to get a decision rule of
the following form. In Einstein notation, it is:
\[
y_t-\bar y = \sum_{i=0}^k\frac{1}{i!}\left[g_{(y^*u)^i}\right]
_{\alpha_1\ldots\alpha_i}
\prod_{j=1}^i\left[\begin{array}{c} y^*_{t-1}-\bar y^*\\ u_t \end{array}\right]
^{\alpha_j}
\]
Note that the ergodic mean will be different from the deterministic
steady state $\bar y$ and thus deviations $y^*_{t-1}-\bar y^*$ will
not be zero in average. This implies that in average we will commit
larger round off errors than if we used the decision rule expressed in
deviations from a point closer to the ergodic mean. Therefore, by
default, Dynare++ recalculates this rule and expresses it in
deviations from the stochastic fix point $y$.
\[
y_t-y = \sum_{i=1}^k\frac{1}{i!}\left[\tilde g_{(y^*u)^i}\right]
_{\alpha_1\ldots\alpha_i}
\prod_{j=1}^i\left[\begin{array}{c} y^*_{t-1}-y^*\\ u_t \end{array}\right]
^{\alpha_j}
\]
Note that since the rule is centralized around its fix point, the
first term (for $i=0$) drops out.
Also note, that this rule mathematically equivalent to the rule
expressed in deviations from the deterministic steady state, and still
it is an approximation about the deterministic steady state. The fact
that it is expressed in deviations from a different point should not
be confused with the algorithm in \ref{multistep_alg}.
This centralization can be avoided by invoking {\tt --no-centralize}
command line option.
\subsection{Taking Steps in Volatility Dimension}
\label{multistep_alg}
For models, where volatility of the exogenous shocks plays a big
role, the approximation about deterministic steady state can be poor,
since the equilibrium dynamics can be very different from the dynamics
in the vicinity of the perfect foresight (deterministic steady state).
Therefore, Dynare++ has on option {\tt --steps} triggering a multistep
algorithm. The algorithm splits the volatility to a given number of
steps. Dynare++ attempts to calculate approximations about fix points
corresponding to these levels of volatility. The problem is that if we
want to calculate higher order approximations about fix points
corresponding to volatilities different from zero (as in the case of
deterministic steady state), then the derivatives of lower orders
depend on derivatives of higher orders with respect to forward looking
variables. The multistep algorithm in each step approximates the
missing higher order derivatives with extrapolations based on the
previous step.
In this way, the approximation of the stochastic fix point and the
derivatives about this fix point are obtained. It is difficult to a
priori decide whether this algorithm yields a better decision
rule. Nothing is guaranteed, and the resulted decision rule should be
checked with a numerical integration. See \ref{checks}.
\subsection{Simulating the Decision Rule}
After some form of a decision rule is calculated, it is simulated to
obtain draws from ergodic (unconditional) distribution of endogenous
variables. The mean and the covariance are reported. There are two
ways how to calculate the mean and the covariance. The first one is to
store all simulated samples and calculate the sample mean and
covariance. The second one is to calculate mean and the covariance in
the real-time not storing the simulated sample. The latter case is
described below (see \ref{rt_simul}).
The stored simulated samples are then used for impulse response
function calculations. For each shock, the realized shocks in these
simulated samples (control simulations) are taken and an impulse is
added and the new realization of shocks is simulated. Then the control
simulation is subtracted from the simulation with the impulse. This is
done for all control simulations and the results are averaged. As the
result, we get an expectation of difference between paths with impulse
and without impulse. In addition, the sample variances are
reported. They might be useful for confidence interval calculations.
For each shock, Dynare++ calculates IRF for two impulses, positive and
negative. Size of an impulse is one standard error of a respective
shock.
The rest of this subsection is divided to three parts giving account
on real-time simulations, conditional simulations, and on the way how
random numbers are generated resp.
\subsubsection{Simulations With Real-Time Statistics}
\label{rt_simul}
When one needs to simulate large samples to get a good estimate of
unconditional mean, simulating the decision rule with statistics
calculated in real-time comes handy. The main reason is that the
storing of all simulated samples may not fit into the available
memory.
The real-time statistics proceed as follows: We model the ergodic
distribution as having normal distribution $y\sim N(\mu,\Sigma)$. Further,
the parameters $\mu$ and $\Sigma$ are modelled as:
\begin{eqnarray*}
\Sigma &\sim& {\rm InvWishart}_\nu(\Lambda)\\
\mu|\Sigma &\sim& N(\bar\mu,\Sigma/\kappa) \\
\end{eqnarray*}
This model of $p(\mu,\Sigma)$ has an advantage of conjugacy, i.e. a
prior distribution has the same form as posterior. This property is
used in the calculation of real-time estimates of $\mu$ and $\Sigma$,
since it suffices to maintain only the parameters of $p(\mu,\Sigma)$
conditional observed draws so far. The parameters are: $\nu$,
$\Lambda$, $\kappa$, and $\bar\mu$.
The mean of $\mu,\Sigma|Y$, where $Y$ are all the draws (simulated
periods) is reported.
\subsubsection{Conditional Distributions}
\label{cond_dist}
Starting with version 1.3.6, Dynare++ calculates variable
distributions $y_t$ conditional on $y_0=\bar y$, where $\bar y$ is the
deterministic steady state. If triggered, Dynare++ simulates a given
number of samples with a given number of periods all starting at
the deterministic steady state. Then for each time $t$, mean
$E[y_t|y_0=\bar y]$ and variances $E[(y_t-E[y_t|y_0=\bar
y])(y_t-E[y_t|y_0=\bar y])^T|y_0=\bar y]$ are reported.
\subsubsection{Random Numbers}
\label{random_numbers}
For generating of the pseudo random numbers, Dynare++ uses Mersenne
twister by Makoto Matsumoto and Takuji Nishimura. Because of the
parallel nature of Dynare++ simulations, each simulated sample gets
its own instance of the twister. Each such instance is seeded before
the simulations are started. This is to prevent additional randomness
implied by the operating system's thread scheduler to interfere with
the pseudo random numbers.
For seeding the individual instances of the Mersenne twister assigned
to each simulated sample the system (C library) random generator is
used. These random generators do not have usually very good
properties, but we use them only to seed the Mersenne twister
instances. The user can set the initial seed of the system random
generator and in this way deterministically choose the seeds of all
instances of the Mersenne twister.
In this way, it is guaranteed that two runs of Dynare++
with the same seed will yield the same results regardless the
operating system's scheduler. The only difference may be caused by a
different round-off errors committed when the same set of samples are
summed in the different order (due to the operating system's scheduler).
\subsection{Numerical Approximation Checks}
\label{checks}
Optionally, Dynare++ can run three kinds of checks for Taylor
approximation errors. All three methods numerically calculate
the residual of the DSGE equations
\[
E[f(g^{**}(g^*(y^*,u),u'),g(y^*,u),y^*,u)|y^*,u]
\]
which must be ideally zero for all $y^*$ and $u$. This integral is
evaluated by either product or Smolyak rule applied to one dimensional
Gauss--Hermite quadrature. The user does not need to care about the
decision. An algorithm yielding higher quadrature level and less
number of evaluations less than a user given maximum is selected.
The three methods differ only by a set of $y^*$ and $u$ where the
residuals are evaluated. These are:
\begin{itemize}
\item The first method calculates the residuals along the shocks for
fixed $y^*$ equal to the fix point. We let all elements of $u$ be
fixed at $0$ but one element, which varies from $-\mu\sigma$ to
$\mu\sigma$, where $\sigma$ is a standard error of the element and
$\mu$ is the user given multiplier. In this way we can see how the
approximation error grows if the fix point is disturbed by a shock of
varying size.
\item The second method calculates the residuals along a simulation
path. A random simulation is run, and at each point the residuals are
reported.
\item The third method calculates the errors on an ellipse of the
state variables $y^*$. The shocks $u$ are always zero. The ellipse is
defined as
\[\{Ax|\; \Vert x\Vert_2=\mu\},\]
where $\mu$ is a user given multiplier, and $AA^T=V$ for $V$ being a
covariance of endogenous variables based on the first order
approximation. The method calculates the residuals at low discrepancy
sequence of points on the ellipse. Both the residuals and the points
are reported.
\end{itemize}
\section{Optimal Policy with Dynare++}
\label{optim}
Starting with version 1.3.2, Dynare++ is able to automatically
generate and then solve the first order conditions for a given
objective and (possibly) forward looking constraints. Since the
constraints can be forward looking, the use of this feature will
mainly be in optimal policy or control.
The only extra thing which needs to be added to the model file is a
specification of the policy's objective. This is done by two keywords,
placed not before parameter settings. If the objective is to maximize
$$E_{t_0}\sum_{t=t_0}^\infty\beta^{t-t_0}\left[\frac{c_t^{1-\eta}}{1-\eta}+
a\frac{g_t^{1-\phi}}{1-\phi}\right],$$
then the keywords will be:
{\small
\begin{verbatim}
planner_objective C^(1-eta)/(1-eta) + a*G^(1-phi)/(1-phi);
planner_discount beta;
\end{verbatim}
}
Dynare++ parses the file and if the two keywords are present, it
automatically derives the first order conditions for the problem. The
first order conditions are put to the form \eqref{focs} and solved. In
this case, the equations in the {\tt model} section are understood as
the constraints (they might come as the first order conditions from
optimizations of other agents) and their number must be less than the
number of endogenous variables.
This section further describes how the optimal policy first order
conditions look like, then discusses some issues with the initial
guess for deterministic steady state, and finally describes how to
simulate Ramsey policy within this framework.
\subsection{First Order Conditions}
Mathematically, the optimization problem looks as follows:
\begin{align}
\max_{\left\{y_\tau\right\}^\infty_t}&E_t
\left[\sum_{\tau=t}^\infty\beta^{\tau-t}b(y_{\tau-1},y_\tau,y_{\tau+1},u_\tau)\right]\notag\\
&\rm{s.t.}\label{planner_optim}\\
&\hskip1cm E^I_\tau\left[f(y_{\tau-1},y_\tau,y_{\tau+1},u_\tau)\right]=0\quad\rm{for\ }
\tau=\ldots,t-1,t,t+1,\ldots\notag
\end{align}
where $E^I$ is an expectation operator over an information set including,
besides all the past, all future realizations of policy's control
variables and distributions of future shocks $u_t\sim
N(0,\Sigma)$. The expectation operator $E$ integrates over an
information including only distributions of $u_t$ (besides the past).
Note that the constraints $f$ take place at all times, and they are
conditioned at the running $\tau$ since the policy knows that the
agents at time $\tau$ will use all the information available at
$\tau$.
The maximization problem can be rewritten using Lagrange multipliers as:
\begin{align}
\max_{y_t}E_t&\left[\sum_{\tau=t}^\infty\beta^{\tau-t}b(y_{\tau-1},y_\tau,y_{\tau+1},u_\tau)\right.\notag\\
&\left.+\sum_{\tau=-\infty}^{\infty}\beta^{\tau-t}\lambda^T_\tau E_\tau^I\left[f(y_{\tau-1},y_\tau,y_{\tau+1},u_\tau)\right]\right],
\label{planner_optim_l}
\end{align}
where $\lambda_t$ is a column vector of Lagrange multipliers.
After some manipulations with compounded expectations over different
information sets, one gets the following first order conditions:
\begin{align}
E_t\left[\vphantom{\frac{\int^(_)}{\int^(\_)}}\right.&\frac{\partial}{\partial y_t}b(y_{t-1},y_t,y_{t+1},u_t)+
\beta L^{+1}\frac{\partial}{\partial y_{t-1}}b(y_{t-1},y_t,y_{t+1},u_t)\notag\\
&+\beta^{-1}\lambda_{t-1}^TL^{-1}\frac{\partial}{\partial y_{t+1}}f(y_{t-1},y_t,y_{t+1},u_t)\notag\\
&+\lambda_t^T\frac{\partial}{\partial y_{t}}f(y_{t-1},y_t,y_{t+1},u_t)\notag\\
&+\beta\lambda_{t+1}^TE_{t+1}\left[L^{+1}\frac{\partial}{\partial y_{t-1}}f(y_{t-1},y_t,y_{t+1},u_t)\right]
\left.\vphantom{\frac{\int^(_)}{\int^(\_)}}\right]
= 0,\label{planner_optim_foc2}
\end{align}
where $L^{+1}$ is one period lead operator, and $L^{-1}$ is one period lag operator.
Dynare++ takes input corresponding to \eqref{planner_optim},
introduces the Lagrange multipliers according to
\eqref{planner_optim_l}, and using its symbolic derivator it compiles
\eqref{planner_optim_foc2}. The system \eqref{planner_optim_foc2} with
the constraints from \eqref{planner_optim_l} is then solved in the
same way as the normal input \eqref{focs}.
\subsection{Initial Guess for Deterministic Steady State}
\label{opt_init}
Solving deterministic steady state of non-linear dynamic systems is
not trivial and the first order conditions for optimal policy add
significant complexity. The {\tt initval} section allows to input the
initial guess of the non-linear solver. It requires that all user
declared endogenous variables be initialized. However, in most cases,
we have no idea what are good initial guesses for the Lagrange
multipliers.
For this reason, Dynare++ calculates an initial guess of Lagrange
multipliers using user provided initial guesses of all other
endogenous variables. It uses the linearity of the Lagrange
multipliers in the \eqref{planner_optim_foc2}. In its static form,
\eqref{planner_optim_foc2} looks as follows:
\begin{align}
&\frac{\partial}{\partial y_t}b(y,y,y,0)+
\beta\frac{\partial}{\partial y_{t-1}}b(y,y,y,0)\notag\\
&+\lambda^T\left[\beta^{-1}\frac{\partial}{\partial y_{t+1}}f(y,y,y,0)
+\frac{\partial}{\partial y_{t}}f(y,y,y,0)
+\beta\frac{\partial}{\partial y_{t-1}}f(y,y,y,0)\right]
= 0\label{planner_optim_static}
\end{align}
The user is required to provide an initial guess of all declared
variables (all $y$). Then \eqref{planner_optim_static} becomes an
overdetermined linear system in $\lambda$, which is solved by means of
the least squares. The closer the initial guess of $y$ is to the exact
solution, the closer are the Lagrange multipliers $\lambda$.
The calculated Lagrange multipliers by the least squares are not used,
if they are set in the {\tt initval} section. In other words, if a
multiplier has been given a value in the {\tt initval} section, then
the value is used, otherwise the calculated value is taken.
For even more difficult problems, Dynare++ generates two Matlab files
calculating a residual of the static system and its derivative. These
can be used in Matlab's {\tt fsolve} or other algorithm to get an
exact solution of the deterministic steady state. See
\ref{output_matlab_scripts} for more details.
Finally, Dynare++ might generate a few auxiliary variables. These are
simple transformations of other variables. They are initialized
automatically and the user usually does not need to care about it.
\subsection{Optimal Ramsey Policy}
\label{ramsey}
Dynare++ solves the optimal policy problem with timeless
perspective. This means that it assumes that the constraints in
\eqref{planner_optim} are valid from the infinite past to infinite
future. Dynare++ calculation of ergodic distribution then assumes that
the policy has been taking place from infinite past.
If some constraints in \eqref{planner_optim} are forward looking, this
will result in some backward looking Lagrange multipliers. Such
multipliers imply possibly time inconsistent policy in the states of
the ``original'' economy, since these backward looking multipliers add
new states to the ``optimized'' economy. In this respect, the timeless
perspective means that there is no fixed initial distribution of such
multipliers, instead, their ergodic distribution is taken.
In contrast, Ramsey optimal policy is started at $t=0$. This means
that the first order conditions at $t=0$ are different than the first
order conditions at $t\geq 1$, which are
\eqref{planner_optim_foc2}. However, it is not difficult to assert
that the first order conditions at $t=0$ are in the form of
\eqref{planner_optim_foc2} if all the backward looking Lagrange
multipliers are set to zeros at period $-1$, i.e. $\lambda_{-1}=0$.
All in all, the solution of \eqref{planner_optim_foc2} calculated by
Dynare++ can be used as a Ramsey optimal policy solution provided that
all the backward looking Lagrange multipliers were set to zeros prior
to the first simulation period. This can be done by setting the
initial state of a simulation path in {\tt dynare\_simul.m}. If this
is applied on the example from \ref{optim_tut}, then we may do the
following in the command prompt:
{\small
\begin{verbatim}
>> load kp1980_2.mat
>> shocks = zeros(1,100);
>> ystart = dyn_ss;
>> ystart(dyn_i_MULT3) = 0;
>> r=dynare_simul('kp1980_2.mat',shocks,ystart);
\end{verbatim}
}
This will simulate the economy if the policy was introduced in the
beginning and no shocks happened.
More information on custom simulations can be obtained by typing:
{\small
\begin{verbatim}
help dynare_simul
\end{verbatim}
}
\section{Running Dynare++}
This section deals with Dynare++ input. The first subsection
\ref{dynpp_opts} provides a list of command line options, next
subsection \ref{dynpp_mod} deals with a format of Dynare++ model file,
and the last subsection discusses incompatibilities between Dynare
Matlab and Dynare++.
\subsection{Command Line Options}
\label{dynpp_opts}
The calling syntax of the Dynare++ is
{\small
\begin{verbatim}
dynare++ [--help] [--version] [options] <model file>
\end{verbatim}
}
\noindent where the model file must be given as the last token and
must include its extension. The model file may include path, in this
case, the path is taken relative to the current directory. Note that
the current directory can be different from the location of {\tt
dynare++} binary.
The options are as follows:
\def\desc#1{\rlap{#1}\kern4cm}
\begin{description}
\item[\desc{\tt --help}] This prints a help message and exits.
\item[\desc{\tt --version}] This prints a version information and
exits.
\item[\desc{\tt --per \it num}] This sets a number of simulated
periods to {\it num} in addition to the burn-in periods. This number
is used when calculating unconditional mean and covariance and for
IRFs. Default is 100.
\item[\desc{\tt --burn \it num}] This sets a number of initial periods
which should be ignored from the statistics. The burn-in periods are
used to eliminate the influence of the starting point when
calculating ergodic distributions or/and impulse response
functions. The number of simulated period given by {\tt --per \it
num} option does not include the number of burn-in
periods. Default is 0.
\item[\desc{\tt --sim \it num}] This sets a number of stochastic
simulations. This number is used when calculating unconditional mean
and covariance and for IRFs. The total sample size for unconditional
mean and covariance is the number of periods times the number of
successful simulations. Note that if a simulation results in {\tt NaN}
or {\tt +-Inf}, then it is thrown away and is not considered for the
mean nor the variance. The same is valid for IRF. Default is 80.
\item[\desc{\tt --rtsim \it num}] This sets a number of stochastic
simulations whose statistics are calculated in the real-time. This
number excludes the burn-in periods set by {\tt --burn \it num}
option. See \ref{rt_simul} for more details. Default is 0, no
simulations.
\item[\desc{\tt --rtper \it num}] This sets a number of simulated
periods per one simulation with real-time statistics to {\it num}. See
\ref{rt_simul} for more details. Default is 0, no simulations.
\item[\desc{\tt --condsim \it num}] This sets a number of stochastic
conditional simulations. See \ref{cond_dist} for more details. Default
is 0, no simulations.
\item[\desc{\tt --condper \it num}] This sets a number of simulated
periods per one conditional simulation. See \ref{cond_dist} for more
details. Default is 0, no simulations.
\item[\desc{\tt --steps \it num}] If the number {\it num} is greater
than 0, this option invokes a multi-step algorithm (see section
\ref{dynpp_calc}), which in the given number of steps calculates fix
points and approximations of the decision rule for increasing
uncertainty. Default is 0, which invokes a standard algorithm for
approximation about deterministic steady state. For more details,
see \ref{multistep_alg}.
\item[\desc{\tt --centralize}] This option causes that the resulting
decision rule is centralized about (in other words: expressed in the
deviations from) the stochastic fix point. The centralized decision
rule is mathematically equivalent but has an advantage of yielding
less numerical errors in average than not centralized decision
rule. By default, the rule is centralized. For more details, see
\ref{dr_form}.
\item[\desc{\tt --no-centralize}] This option causes that the
resulting decision rule is not centralized about (in other words:
expressed in the deviations from) the stochastic fix point. By
default, the rule is centralized. For more details, see
\ref{dr_form}.
This option has no effect if the number of steps given by {\tt
--steps} is greater than 0. In this case, the rule is always
centralized.
\item[\desc{\tt --prefix \it string}] This sets a common prefix of
variables in the output MAT file. Default is {\tt dyn}.
\item[\desc{\tt --seed \it num}] This sets an initial seed for the
random generator providing seed to generators for each sample. See
\ref{random_numbers} for more details. Default is 934098.
\item[\desc{\tt --order \it num}] This sets the order of approximation
and overrides the {\tt order} statement in the model file. There is no
default.
\item[\desc{\tt --threads \it num}] This sets a number of parallel
threads. Complex evaluations of Faa Di Bruno formulas, simulations and
numerical integration can be parallelized, Dynare++ exploits this
advantage. You have to have a hardware support for this, otherwise
there is no gain from the parallelization. The default value is the number of
logical processors present on the machine, divided by 2.
\item[\desc{\tt --ss-tol \it float}] This sets the tolerance of the
non-linear solver of deterministic steady state to {\it float}. It is
in $\Vert\cdot\Vert_\infty$ norm, i.e. the algorithm is considered as
converged when a maximum absolute residual is less than the
tolerance. Default is $10^{-13}$.
\item[\desc{\tt --check \it pPeEsS}] This selects types of residual
checking to be performed. See section \ref{checks} for details. The
string consisting of the letters ``pPeEsS'' governs the selection. The
upper-case letters switch a check on, the lower-case letters
off. ``P'' stands for checking along a simulation path, ``E'' stands
for checking on ellipse, and finally ``S'' stands for checking along
the shocks. It is possible to choose more than one type of check. The
default behavior is that no checking is performed.
\item[\desc{\tt --check-evals \it num}] This sets a maximum number of
evaluations per one re\-sidual. The actual value depends on the selected
algorithm for the integral evaluation. The algorithm can be either
product or Smolyak quadrature and is chosen so that the actual number
of evaluations would be minimal with maximal level of
quadrature. Default is 1000.
\item[\desc{\tt --check-num \it num}] This sets a number of checked
points in a residual check. One input value $num$ is used for all
three types of checks in the following way:
\begin{itemize}
\item For checks along the simulation, the number of simulated periods
is $10\cdot num$
\item For checks on ellipse, the number of points on ellipse is $10\cdot num$
\item For checks along the shocks, the number of checked points
corresponding to shocks from $0$ to $\mu\sigma$ (see \ref{checks}) is
$num$.
\end{itemize}
Default is 10.
\item[\desc{\tt --check-scale \it float}] This sets the scaling factor
$\mu$ for checking on ellipse to $0.5\cdot float$ and scaling factor
$\mu$ for checking along shocks to $float$. See section
\ref{checks}. Default is 2.0.
\item[\desc{\tt --no-irfs}] This suppresses IRF calculations. Default
is to calculate IRFs for all shocks.
\item[\desc{\tt --irfs}] This triggers IRF calculations. If there are
no shock names following the {\tt --irfs} option, then IRFs for all
shocks are calculated, otherwise see below. Default is to calculate
IRFs for all shocks.
\item[\desc{\tt --irfs \it shocklist}] This triggers IRF calculations
only for the listed shocks. The {\it shocklist} is a space separated
list of exogenous variables for which the IRFs will be
calculated. Default is to calculate IRFs for all shocks.
\end{description}
The following are a few examples:
{\small
\begin{verbatim}
dynare++ --sim 300 --per 50 blah.mod
dynare++ --check PE --check-num 15 --check-evals 500 blah.dyn
dynare++ --steps 5 --check S --check-scale 3 blahblah.mod
\end{verbatim}
}
The first one sets the number of periods for IRF to 50, and sets a sample
size for unconditional mean and covariance calculations to 6000. The
second one checks the decision rule along a simulation path having 150
periods and on ellipse at 150 points performing at most 500 evaluations
per one residual. The third one solves the model in five steps and
checks the rule along all the shocks from $-3\sigma$ to $3\sigma$ in
$2*10+1$ steps (10 for negative, 10 for positive and 1 for at zero).
\subsection{Dynare++ Model File}
\label{dynpp_mod}
In its strictest form, Dynare++ solves the following mathematical problem:
\begin{equation}\label{basic_form}
E_t[f(y^{**}_{t+1},y_t,y^*_{t-1},u_t)]=0
\end{equation}
This problem is input either directly, or it is an output of Dynare++
routines calculating first order conditions of the optimal policy
problem. In either case, Dynare++ performs necessary and
mathematically correct substitutions to put the user specified problem
to the \eqref{basic_form} form, which goes to Dynare++ solver. The
following discusses a few timing issues:
\begin{itemize}
\item Endogenous variables can occur, starting from version 1.3.4, at
times after $t+1$. If so, an equation containing such occurrence is
broken to non-linear parts, and new equations and new auxiliary
variables are automatically generated only for the non-linear terms
containing the occurrence. Note that shifting such terms to time $t+1$
may add occurrences of some other variables (involved in the terms) at
times before $t-1$ implying addition of auxiliary variables to bring
those variables to $t-1$.
\item Variables declared as shocks may occur also at arbitrary
times. If before $t$, additional endogenous variables are used to
bring them to time $t$. If after $t$, then similar method is used as
for endogenous variables occurring after $t+1$.
\item There is no constraint on variables occurring at both times
$t+1$ (or later) and $t-1$ (or earlier). Virtually, all variables can
occur at arbitrary times.
\item Endogenous variables can occur at times before $t-1$. If so,
additional endogenous variables are added for all lags between the
variable and $t-1$.
\item Dynare++ applies the operator $E_t$ to all occurrences at time
$t+1$. The realization of $u_t$ is included in the information set of
$E_t$. See an explanation of Dynare++ timing on page \pageref{timing}.
\end{itemize}
The model equations are formulated in the same way as in Matlab
Dynare. The time indexes different from $t$ are put to round
parenthesis in this way: {\tt C(-1)}, {\tt C}, {\tt C(+1)}.
The mathematical expressions can use the following functions and operators:
\begin{itemize}
\item binary {\tt + - * / \verb|^|}
\item unary plus and minus minus as in {\tt a = -3;} and {\tt a = +3;} resp.
\item unary mathematical functions: {\tt log exp sin cos tan
sqrt}, whe\-re the logarithm has a natural base
\item symbolic differentiation operator {\tt diff(expr,symbol)}, where
{\tt expr} is a mathematical expression and {\tt symbol} is a unary
symbol (a variable or a parameter); for example {\tt
diff(A*K(-1)\verb|^|alpha*L\verb|^|(1-alpha),K(-1))} is internally expanded as
{\tt A*alpha*K(-1)\verb|^|(alpha-1)*L\verb|^|(1-alpha)}
\item unary error function and complementary error function: {\tt erf}
and {\tt erfc} defined as
\begin{eqnarray*}
erf(x) &= \frac{2}{\sqrt{\pi}}\int_0^x e^{-t^2}{\rm d}t\\
erfc(x)&= \frac{2}{\sqrt{\pi}}\int_x^\infty e^{-t^2}{\rm d}t
\end{eqnarray*}
\end{itemize}
The model file can contain user comments. Their usage can be
understood from the following piece of the model file:
{\small
\begin{verbatim}
P*C^(-gamma) = // line continues until semicolon
beta*C(+1)^(-gamma)*(P(+1)+Y(+1)); // asset price
// choose dividend process: (un)comment what you want
Y/Y_SS = (Y(-1)/Y_SS)^rho*exp(EPS);
/*
Y-Y_SS = rho*(Y(-1)-Y_SS)+EPS;
*/
\end{verbatim}
}
\subsection{Incompatibilities with Matlab Dynare}
This section provides a list of incompatibilities between a model file
for Dy\-na\-re++ and Matlab Dynare. These must be considered when a model
file for Matlab Dynare is being migrated to Dynare++. The list is the
following:
\begin{itemize}
\item There is no {\tt periods} keyword.
\item The parameters cannot be lagged or leaded, I think that Dynare
Matlab allows it, but the semantics is the same (parameter is a
constant).
\item There are no commands like {\tt steady}, {\tt check}, {\tt
simul}, {\tt stoch\_simul}, etc.
\item There are no sections like {\tt estimated\_params}, {\tt
var\_obs}, etc.
\item The variance-covariance matrix of endogenous shocks is given by
{\tt vcov} matrix in Dynare++. An example follows. Starting from
version 1.3.5, it is possible for vcov to be positive semi-definite
matrix.
{\small
\begin{verbatim}
vcov = [
0.05 0 0 0;
0 0.025 0 0;
0 0 0.05 0;
0 0 0 0.025
];
\end{verbatim}
}
\end{itemize}
\section{Dynare++ Output}
There are three output files; a data file in MAT-4 format containing
the output data (\ref{matfile}), a journal text file containing an
information about the Dynare++ run (\ref{journalfile}), and a dump
file (\ref{dumpfile}). Further, Dynare++ generates two Matlab script
files, which calculate a residual and the first derivative of the
residual of the static system (\ref{output_matlab_scripts}). These are
useful when calculating the deterministic steady state outside
Dynare++.
Note that all output files are created in the current directory of
the Dynare++ process. This can be different from the directory where
the Dynare++ binary is located and different from the directory where
the model file is located.
Before all, we need to understand what variables are automatically
generated in Dynare++.
\subsection{Auxiliary Variables}
\label{aux_var}
Besides the endogenous variables declared in {\tt var} section,
Dynare++ might automatically add the following endogenous variables:
\halign{\vrule width0pt height14pt{\tt #}\hfil & \kern 3mm%
\vtop{\rightskip=0pt plus 5mm\noindent\hsize=9.5cm #}\cr
MULT{\it n}& A Lagrange multiplier of the optimal policy problem
associated with a constraint number {\it n} starting from zero.\cr
AUX\_{\it n1}\_{\it n2}\_{\it n3}& An auxiliary variable associated
with the last term in equation \eqref{planner_optim_foc2}. Since the
term is under $E_{t+k}$, we need the auxiliary variable be put back
in time. {\it n1} is a variable number starting from 0 in the declared
order with respect to which the term was differentiated, {\it n2} is a
number of constraint starting from 0, and finally {\it n3} is $k$
(time shift of the term).\cr
{\it endovar}\_p{\it K}& An auxiliary variable for bringing an
endogenous variable {\it endovar} back in time by $K$ periods. The
semantics of this variables is {\tt {\it endovar}\_p{\it K} = {\it
endovar}(+{\it K})}.\cr
{\it endovar}\_m{\it K}& An auxiliary variable for bringing an
endogenous variable {\it endovar} forward in time by $K$ periods. The
semantics of this variables is {\tt {\it endovar}\_m{\it K} = {\it
endovar}(-{\it K})}.\cr
{\it exovar}\_e& An auxiliary endogenous variable made equal to the
exogenous variable to allow for a semantical occurrence of the
exogenous variable at time other than $t$. The semantics of this
variables is {\tt {\it exovar}\_e = {\it exovar}}.\cr
AUXLD\_{\it n1}\_{\it n2}\_{\it n3}& An auxiliary variable for
bringing a non-linear term containing an occurrence of a variable
after $t+1$ to time $t+1$. {\it n1} is an equation number starting
from 0, {\it n2} is the non-linear sub-term number in the equation
starting from 0. {\it n3} is a time shift. For example, if the first
equation is the following:
\begin{verbatim}
X - Y*W(+1) + W(+2)*Z(+4) = 0;
\end{verbatim}
then it will be expanded as:
\begin{verbatim}
X - Y*W(+1) + AUXLD_0_2_3(+1) = 0;
AUXLD_0_2_1 = W(-1)*Z(+1);
AUXLD_0_2_2 = AUXLD_0_2_1(+1);
AUXLD_0_2_3 = AUXLD_0_2_2(+1);
\end{verbatim}
\cr
}
\subsection{MAT File}
\label{matfile}
The contents of the data file is depicted below. We
assume that the prefix is {\tt dyn}.
\halign{\vrule width0pt height14pt{\tt #}\hfil & \kern 3mm%
\vtop{\rightskip=0pt plus 5mm\noindent\hsize=7.5cm #}\cr
dyn\_nstat& Scalar. A number of static variables
(those occurring only at time $t$).\cr
dyn\_npred & Scalar. A number of variables occurring
at time $t-1$ and not at $t+1$.\cr
dyn\_nboth & Scalar. A number of variables occurring
at $t+1$ and $t-1$.\cr
dyn\_nforw & Scalar. A number of variables occurring
at $t+1$ and not at $t-1$.\cr
dyn\_vars & Column vector of endogenous variable
names in Dy\-na\-re++ internal ordering.\cr
dyn\_i\_{\it endovar} & Scalar. Index of a variable
named {\it endovar} in the {\tt dyn\_vars}.\cr
dyn\_shocks & Column vector of exogenous variable
names.\cr
dyn\_i\_{\it exovar} & Scalar. Index of a shock
named {\it exovar} in the {\tt dyn\_shocks}.\cr
dyn\_state\_vars & Column vector of state variables,
these are stacked variables counted by {\tt dyn\_\-npred}, {\tt
dyn\_\-nboth} and shocks.\cr
dyn\_vcov\_exo & Matrix $nexo\times nexo$. The
variance-covariance matrix of exogenous shocks as input in the model
file. The ordering is given by {\tt dyn\_shocks}.\cr
dyn\_mean & Column vector $nendo\times 1$. The
unconditional mean of endogenous variables. The ordering is given by
{\tt dyn\_vars}.\cr
dyn\_vcov & Matrix $nendo\times nendo$. The
unconditional covariance of endogenous variables. The ordering is given
by {\tt dyn\_vars}.\cr
dyn\_rt\_mean & Column vector $nendo\times 1$. The unconditional mean
of endogenous variables estimated in real-time. See
\ref{rt_simul}. The ordering is given by {\tt dyn\_vars}.\cr
dyn\_rt\_vcov & Matrix $nendo\times nendo$. The unconditional
covariance of endogenous variables estimated in real-time. See \ref{rt_simul}. The
ordering is given by {\tt dyn\_vars}.\cr
dyn\_cond\_mean & Matrix $nendo\times nper$. The rows correspond to
endogenous variables in the ordering of {\tt dyn\_vars}, the columns
to periods. If $t$ is a period (starting with 1), then $t$-th column
is $E[y_t|y_0=\bar y]$. See \ref{cond_dist}.\cr
dyn\_cond\_variance & Matrix $nendo\times nper$. The rows correspond
to endogenous variables in the ordering of {\tt dyn\_vars}, the
columns to periods. If $t$ is a period (starting with 1), then $t$-th
column are the variances of $y_t|y_0=\bar y$. See \ref{cond_dist}.\cr
dyn\_ss & Column vector $nendo\times 1$. The fix
point of the resulting approximation of the decision rule.\cr
dyn\_g\_{\it order} & Matrix $nendo\times ?$. A
derivative of the decision rule of the {\it order} multiplied by
$1/order!$. The rows correspond to endogenous variables in the
ordering of {\tt dyn\_vars}. The columns correspond to a
multidimensional index going through {\tt dyn\_state\_vars}. The data
is folded (all symmetrical derivatives are stored only once).\cr
dyn\_steady\_states & Matrix $nendo\times
nsteps+1$. A list of fix points at which the multi-step algorithm
calculated approximations. The rows correspond to endogenous variables
and are ordered by {\tt dyn\_vars}, the columns correspond to the
steps. The first column is always the deterministic steady state.\cr
dyn\_irfp\_{\it exovar}\_mean & Matrix
$nendo\times nper$. Positive impulse response to a shock named {\it
exovar}. The row ordering is given by {\tt dyn\_vars}. The columns
correspond to periods.\cr
dyn\_irfp\_{\it exovar}\_var & Matrix
$nendo\times nper$. The variances of positive impulse response
functions.\cr
dyn\_irfm\_{\it exovar}\_mean & Same as {\tt
dyn\_irfp\_}{\it exovar}{\tt \_mean} but for negative impulse.\cr
dyn\_irfp\_{\it exovar}\_var & Same as {\tt
dyn\_irfp\_}{\it exovar}{\tt \_var} but for negative impulse.\cr
dyn\_simul\_points & A simulation path along which the check was
done. Rows correspond to endogenous variables, columns to
periods. Appears only if {\tt --check P}.\cr
dyn\_simul\_errors & Errors along {\tt
dyn\_simul\_points}. The rows correspond to equations as stated in the
model file, the columns to the periods. Appears only if {\tt --check
P}\cr
dyn\_ellipse\_points & A set of points on the ellipse at which the
approximation was checked. Rows correspond to state endogenous
variables (the upper part of {\tt dyn\_state\_vars}, this means
without shocks), and columns correspond to periods. Appears only if
{\tt --check E}.\cr
dyn\_ellipse\_errors & Errors on the ellipse points {\tt
dyn\_ellipse\_points}. The rows correspond to the equations as stated
in the model file, columns to periods. Appears only if {\tt --check
E}.\cr
dyn\_shock\_{\it exovar}\_errors& Errors along a shock named {\it
exovar}. The rows correspond to the equations as stated in the model
file. There are $2m+1$ columns, the middle column is the error at zero
shock. The columns to the left correspond to negative values, columns
to the right to positive. Appears only if {\tt --check S}.\cr
}
\subsection{Journal File}
\label{journalfile}
The journal file provides information on resources usage during the
run and gives some informative messages. The journal file is a text
file, it is organized in single line records. The format of records is
documented in a header of the journal file.
The journal file should be consulted in the following circumstances:
\begin{itemize}
\item Something goes wrong. For example, if a model is not
Blanchard--Kahn stable, then the eigenvalues are dumped to the journal
file.
If the unconditional covariance matrix {\tt dyn\_vcov} is NaN, then
from the journal file you will know that all the simulations had to be
thrown away due to occurrence of NaN or Inf. This is caused by
non-stationarity of the resulting decision rule.
If Dynare++ crashes, the journal file can be helpful for guessing a
point where it crashed.
\item You are impatient. You might be looking at the journal file
during the run in order to have a better estimate about the time when
the calculations are finished. In Unix, I use a command {\tt tail -f
blah.jnl}.\footnote{This helps to develop one of the three
programmer's virtues: {\it impatience}. The other two are {\it
laziness} and {\it hubris}; according to Larry Wall.}
\item Heavy swapping. If the physical memory is not
sufficient, an operating system starts swapping memory pages with a
disk. If this is the case, the journal file can be consulted for
information on memory consumption and swapping activity.
\item Not sure what Dynare++ is doing. If so, read the journal file,
which contains a detailed record on what was calculated, simulated
etc.
\end{itemize}
\subsection{Dump File}
\label{dumpfile}
The dump file is always created with the suffix {\tt .dump}. It is a
text file which takes a form of a model file. It sets the parameter
values which were used, it has the initval section setting the values
which were finally used, and mainly it has a model section of all
equations with all substitutions and formed the first order conditions
of the planner.
The dump file serves for debugging purposes, since it contains the
mathematical problem which is being solved by dynare++.
\subsection{Matlab Scripts for Steady State Calculations}
\label{output_matlab_scripts}
This section describes two Matlab scripts, which are useful when
calculating the deterministic steady state outside Dynare++. The
scripts are created by Dynare++ as soon as an input file is parsed,
that is before any calculations.
The first Matlab script having a name {\tt {\it modname}\_f.m} for
given parameters values and given all endogenous variables $y$
calculates a residual of the static system. Supposing the model is in
the form of \eqref{focs}, the script calculates a vector:
\[
f(y,y,y,0)
\]
The second script having a name {\tt {\it modname}\_ff.m} calculates a matrix:
\[
\frac{\partial}{\partial y}f(y,y,y,0)
\]
Both scripts take two arguments. The first is a vector of parameter
values ordered in the same ordering as declared in the model file. The
second is a vector of all endogenous variables at which the evaluation
is performed. These endogenous variables also include auxiliary
variables automatically added by Dynare++ and Lagrange multipliers if
an optimal policy problem is solved. If no endogenous variable has not
been added by Dynare++, then the ordering is the same as the ordering
in declaration in the model file. If some endogenous variables have
been added, then the ordering can be read from comments close to the
top of either two files.
For example, if we want to calculate the deterministic steady state of
the {\tt kp1980.dyn} model, we need to do the following:
\begin{enumerate}
\item Run Dynare++ with {\tt kp1980.dyn}, no matter if the calculation
has not been finished, important output are the two Matlab scripts
created just in the beginning.
\item Consult file {\tt kp1980\_f.m}\ to get the ordering of parameters
and all endogenous variables.
\item Create a vector {\tt p} with the parameter values in the ordering
\item Create a vector {\tt init\_y} with the initial guess for the
Matlab solver {\tt fsolve}
\item Create a simple Matlab function called {\tt kp1980\_fsolve.m}\
returning the residual and Jacobian:
{\small
\begin{verbatim}
function [r, J] = kp1980_fsolve(p, y)
r = kp1980_f(p, y);
J = kp1980_ff(p, y);
\end{verbatim}
}
\item In the Matlab prompt, run the following:
{\small
\begin{verbatim}
opt=optimset('Jacobian','on','Display','iter');
y=fsolve(@(y) kp1980_fsolve(p,y), init_y, opt);
\end{verbatim}
}
\end{enumerate}
\subsection{Custom Simulations}
\label{custom}
When Dynare++ run is finished it dumps the derivatives of the
calculated decision rule to the MAT file. The derivatives can be used
for a construction of the decision rule and custom simulations can be
run. This is done by {\tt dynare\_simul.m} M-file in Matlab. It reads
the derivatives and simulates the decision rule with provided
realization of shocks.
All the necessary documentation can be viewed by the command:
{\small
\begin{verbatim}
help dynare_simul
\end{verbatim}
}
\end{document}
dynare++/integ/Makefile.am deleted 100644 → 0
View file @ 53ae5493
SUBDIRS = cc src testing
dynare++/integ/cc/Makefile.am deleted 100644 → 0
View file @ 53ae5493
noinst_LIBRARIES = libinteg.a
libinteg_a_SOURCES = \
quadrature.cc \
quadrature.hh \
quasi_mcarlo.cc \
quasi_mcarlo.hh \
product.cc \
product.hh \
smolyak.cc \
smolyak.hh \
vector_function.cc \
vector_function.hh \
precalc_quadrature.hh
libinteg_a_CPPFLAGS = -I../../sylv/cc -I../../utils/cc -I../../tl/cc -I$(top_srcdir)/mex/sources
libinteg_a_CXXFLAGS = $(AM_CXXFLAGS) $(THREAD_CXXFLAGS)
dynare++/integ/cc/precalc_quadrature.hh deleted 100644 → 0
View file @ 53ae5493
/*
* Copyright © 2005 Ondra Kamenik
* Copyright © 2019 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/>.
*/
// The file contains one dimensional quadrature points and weights for a few
// quadratures. The format of data is clear. There is a class
// OneDPrecalcQuadrature which implements an interface OneDQuadrature using the
// data of this format.
// Gauss-Hermite quadrature; prefix gh
// Number of levels
static const int gh_num_levels = 26;
// Number of points in each level
static const int gh_num_points[] = {
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
30, 32, 40, 50, 60, 64
};
// Weights, starting with the first level
static const double gh_weights[] = {
// weights 1 = √π
1.77245385090551588191942755656782537698745727539062,
// weights 2
0.886226925452758013649083741671e+00,
0.886226925452758013649083741671e+00,
// weights 3
0.295408975150919337883027913890e+00,
0.118163590060367735153211165556e+01,
0.295408975150919337883027913890e+00,
// weights 4
0.813128354472451771430345571899e-01,
0.804914090005512836506049184481e+00,
0.804914090005512836506049184481e+00,
0.813128354472451771430345571899e-01,
// weights 5
0.199532420590459132077434585942e-01,
0.393619323152241159828495620852e+00,
0.945308720482941881225689324449e+00,
0.393619323152241159828495620852e+00,
0.199532420590459132077434585942e-01,
// weights 6
0.453000990550884564085747256463e-02,
0.157067320322856643916311563508e+00,
0.724629595224392524091914705598e+00,
0.724629595224392524091914705598e+00,
0.157067320322856643916311563508e+00,
0.453000990550884564085747256463e-02,
// weights 7
0.971781245099519154149424255939e-03,
0.545155828191270305921785688417e-01,
0.425607252610127800520317466666e+00,
0.810264617556807326764876563813e+00,
0.425607252610127800520317466666e+00,
0.545155828191270305921785688417e-01,
0.971781245099519154149424255939e-03,
// weights 8
0.199604072211367619206090452544e-03,
0.170779830074134754562030564364e-01,
0.207802325814891879543258620286e+00,
0.661147012558241291030415974496e+00,
0.661147012558241291030415974496e+00,
0.207802325814891879543258620286e+00,
0.170779830074134754562030564364e-01,
0.199604072211367619206090452544e-03,
// weights 9
0.396069772632643819045862946425e-04,
0.494362427553694721722456597763e-02,
0.884745273943765732879751147476e-01,
0.432651559002555750199812112956e+00,
0.720235215606050957124334723389e+00,
0.432651559002555750199812112956e+00,
0.884745273943765732879751147476e-01,
0.494362427553694721722456597763e-02,
0.396069772632643819045862946425e-04,
// weights 10
0.764043285523262062915936785960e-05,
0.134364574678123269220156558585e-02,
0.338743944554810631361647312776e-01,
0.240138611082314686416523295006e+00,
0.610862633735325798783564990433e+00,
0.610862633735325798783564990433e+00,
0.240138611082314686416523295006e+00,
0.338743944554810631361647312776e-01,
0.134364574678123269220156558585e-02,
0.764043285523262062915936785960e-05,
// weights 11
0.143956039371425822033088366032e-05,
0.346819466323345510643413772940e-03,
0.119113954449115324503874202916e-01,
0.117227875167708503381788649308e+00,
0.429359752356125028446073598601e+00,
0.654759286914591779203940657627e+00,
0.429359752356125028446073598601e+00,
0.117227875167708503381788649308e+00,
0.119113954449115324503874202916e-01,
0.346819466323345510643413772940e-03,
0.143956039371425822033088366032e-05,
// weights 12
0.265855168435630160602311400877e-06,
0.857368704358785865456906323153e-04,
0.390539058462906185999438432620e-02,
0.516079856158839299918734423606e-01,
0.260492310264161129233396139765e+00,
0.570135236262479578347113482275e+00,
0.570135236262479578347113482275e+00,
0.260492310264161129233396139765e+00,
0.516079856158839299918734423606e-01,
0.390539058462906185999438432620e-02,
0.857368704358785865456906323153e-04,
0.265855168435630160602311400877e-06,
// weights 13
0.482573185007313108834997332342e-07,
0.204303604027070731248669432937e-04,
0.120745999271938594730924899224e-02,
0.208627752961699392166033805050e-01,
0.140323320687023437762792268873e+00,
0.421616296898543221746893558568e+00,
0.604393187921161642342099068579e+00,
0.421616296898543221746893558568e+00,
0.140323320687023437762792268873e+00,
0.208627752961699392166033805050e-01,
0.120745999271938594730924899224e-02,
0.204303604027070731248669432937e-04,
0.482573185007313108834997332342e-07,
// weights 14
0.862859116812515794532041783429e-08,
0.471648435501891674887688950105e-05,
0.355092613551923610483661076691e-03,
0.785005472645794431048644334608e-02,
0.685055342234652055387163312367e-01,
0.273105609064246603352569187026e+00,
0.536405909712090149794921296776e+00,
0.536405909712090149794921296776e+00,
0.273105609064246603352569187026e+00,
0.685055342234652055387163312367e-01,
0.785005472645794431048644334608e-02,
0.355092613551923610483661076691e-03,
0.471648435501891674887688950105e-05,
0.862859116812515794532041783429e-08,
// weights 15
0.152247580425351702016062666965e-08,
0.105911554771106663577520791055e-05,
0.100004441232499868127296736177e-03,
0.277806884291277589607887049229e-02,
0.307800338725460822286814158758e-01,
0.158488915795935746883839384960e+00,
0.412028687498898627025891079568e+00,
0.564100308726417532852625797340e+00,
0.412028687498898627025891079568e+00,
0.158488915795935746883839384960e+00,
0.307800338725460822286814158758e-01,
0.277806884291277589607887049229e-02,
0.100004441232499868127296736177e-03,
0.105911554771106663577520791055e-05,
0.152247580425351702016062666965e-08,
// weights 16
0.265480747401118224470926366050e-09,
0.232098084486521065338749423185e-06,
0.271186009253788151201891432244e-04,
0.932284008624180529914277305537e-03,
0.128803115355099736834642999312e-01,
0.838100413989858294154207349001e-01,
0.280647458528533675369463335380e+00,
0.507929479016613741913517341791e+00,
0.507929479016613741913517341791e+00,
0.280647458528533675369463335380e+00,
0.838100413989858294154207349001e-01,
0.128803115355099736834642999312e-01,
0.932284008624180529914277305537e-03,
0.271186009253788151201891432244e-04,
0.232098084486521065338749423185e-06,
0.265480747401118224470926366050e-09,
// weights 17
0.458057893079863330580889281222e-10,
0.497707898163079405227863353715e-07,
0.711228914002130958353327376218e-05,
0.298643286697753041151336643059e-03,
0.506734995762753791170069495879e-02,
0.409200341495762798094994877854e-01,
0.172648297670097079217645196219e+00,
0.401826469470411956577635085257e+00,
0.530917937624863560331883103379e+00,
0.401826469470411956577635085257e+00,
0.172648297670097079217645196219e+00,
0.409200341495762798094994877854e-01,
0.506734995762753791170069495879e-02,
0.298643286697753041151336643059e-03,
0.711228914002130958353327376218e-05,
0.497707898163079405227863353715e-07,
0.458057893079863330580889281222e-10,
// weights 18
0.782819977211589102925147471012e-11,
0.104672057957920824443559608435e-07,
0.181065448109343040959702385911e-05,
0.918112686792940352914675407371e-04,
0.188852263026841789438175325426e-02,
0.186400423875446519219315221973e-01,
0.973017476413154293308537234155e-01,
0.284807285669979578595606820713e+00,
0.483495694725455552876410522141e+00,
0.483495694725455552876410522141e+00,
0.284807285669979578595606820713e+00,
0.973017476413154293308537234155e-01,
0.186400423875446519219315221973e-01,
0.188852263026841789438175325426e-02,
0.918112686792940352914675407371e-04,
0.181065448109343040959702385911e-05,
0.104672057957920824443559608435e-07,
0.782819977211589102925147471012e-11,
// weights 19
0.132629709449851575185289154385e-11,
0.216305100986355475019693077221e-08,
0.448824314722312295179447915594e-06,
0.272091977631616257711941025214e-04,
0.670877521407181106194696282100e-03,
0.798886677772299020922211491861e-02,
0.508103869090520673569908110358e-01,
0.183632701306997074156148485766e+00,
0.391608988613030244504042313621e+00,
0.502974888276186530840731361096e+00,
0.391608988613030244504042313621e+00,
0.183632701306997074156148485766e+00,
0.508103869090520673569908110358e-01,
0.798886677772299020922211491861e-02,
0.670877521407181106194696282100e-03,
0.272091977631616257711941025214e-04,
0.448824314722312295179447915594e-06,
0.216305100986355475019693077221e-08,
0.132629709449851575185289154385e-11,
// weights 20
0.222939364553415129252250061603e-12,
0.439934099227318055362885145547e-09,
0.108606937076928169399952456345e-06,
0.780255647853206369414599199965e-05,
0.228338636016353967257145917963e-03,
0.324377334223786183218324713235e-02,
0.248105208874636108821649525589e-01,
0.109017206020023320013755033535e+00,
0.286675505362834129719659706228e+00,
0.462243669600610089650328639861e+00,
0.462243669600610089650328639861e+00,
0.286675505362834129719659706228e+00,
0.109017206020023320013755033535e+00,
0.248105208874636108821649525589e-01,
0.324377334223786183218324713235e-02,
0.228338636016353967257145917963e-03,
0.780255647853206369414599199965e-05,
0.108606937076928169399952456345e-06,
0.439934099227318055362885145547e-09,
0.222939364553415129252250061603e-12,
// weights 30
0.290825470013122622941102747365e-20,
0.281033360275090370876277491534e-16,
0.287860708054870606219239791142e-13,
0.810618629746304420399344796173e-11,
0.917858042437852820850075742492e-09,
0.510852245077594627738963204403e-07,
0.157909488732471028834638794022e-05,
0.293872522892298764150118423412e-04,
0.348310124318685523420995323183e-03,
0.273792247306765846298942568953e-02,
0.147038297048266835152773557787e-01,
0.551441768702342511680754948183e-01,
0.146735847540890099751693643152e+00,
0.280130930839212667413493211293e+00,
0.386394889541813862555601849165e+00,
0.386394889541813862555601849165e+00,
0.280130930839212667413493211293e+00,
0.146735847540890099751693643152e+00,
0.551441768702342511680754948183e-01,
0.147038297048266835152773557787e-01,
0.273792247306765846298942568953e-02,
0.348310124318685523420995323183e-03,
0.293872522892298764150118423412e-04,
0.157909488732471028834638794022e-05,
0.510852245077594627738963204403e-07,
0.917858042437852820850075742492e-09,
0.810618629746304420399344796173e-11,
0.287860708054870606219239791142e-13,
0.281033360275090370876277491534e-16,
0.290825470013122622941102747365e-20,
// weights 32
0.731067642736e-22,
0.923173653649e-18,
0.119734401709e-14,
0.421501021125e-12,
0.593329146300e-10,
0.409883216476e-08,
0.157416779254e-06,
0.365058512955e-05,
0.541658406172e-04,
0.536268365526e-03,
0.365489032664e-02,
0.175534288315e-01,
0.604581309557e-01,
0.151269734076e+00,
0.277458142302e+00,
0.375238352592e+00,
0.375238352592e+00,
0.277458142302e+00,
0.151269734076e+00,
0.604581309557e-01,
0.175534288315e-01,
0.365489032664e-02,
0.536268365526e-03,
0.541658406172e-04,
0.365058512955e-05,
0.157416779254e-06,
0.409883216476e-08,
0.593329146300e-10,
0.421501021125e-12,
0.119734401709e-14,
0.923173653649e-18,
0.731067642736e-22,
// weights 40
0.259104371384e-28,
0.854405696375e-24,
0.256759336540e-20,
0.198918101211e-17,
0.600835878947e-15,
0.880570764518e-13,
0.715652805267e-11,
0.352562079135e-09,
0.112123608322e-07,
0.241114416359e-06,
0.363157615067e-05,
0.393693398108e-04,
0.313853594540e-03,
0.187149682959e-02,
0.846088800823e-02,
0.293125655361e-01,
0.784746058652e-01,
0.163378732713e+00,
0.265728251876e+00,
0.338643277425e+00,
0.338643277425e+00,
0.265728251876e+00,
0.163378732713e+00,
0.784746058652e-01,
0.293125655361e-01,
0.846088800823e-02,
0.187149682959e-02,
0.313853594540e-03,
0.393693398108e-04,
0.363157615067e-05,
0.241114416359e-06,
0.112123608322e-07,
0.352562079135e-09,
0.715652805267e-11,
0.880570764518e-13,
0.600835878947e-15,
0.198918101211e-17,
0.256759336540e-20,
0.854405696375e-24,
0.259104371384e-28,
// weights 50
0.183379404857e-36,
0.167380166790e-31,
0.121524412340e-27,
0.213765830835e-24,
0.141709359957e-21,
0.447098436530e-19,
0.774238295702e-17,
0.809426189344e-15,
0.546594403180e-13,
0.250665552389e-11,
0.811187736448e-10,
0.190904054379e-08,
0.334679340401e-07,
0.445702996680e-06,
0.458168270794e-05,
0.368401905377e-04,
0.234269892109e-03,
0.118901178175e-02,
0.485326382616e-02,
0.160319410684e-01,
0.430791591566e-01,
0.945489354768e-01,
0.170032455676e+00,
0.251130856331e+00,
0.305085129203e+00,
0.305085129203e+00,
0.251130856331e+00,
0.170032455676e+00,
0.945489354768e-01,
0.430791591566e-01,
0.160319410684e-01,
0.485326382616e-02,
0.118901178175e-02,
0.234269892109e-03,
0.368401905377e-04,
0.458168270794e-05,
0.445702996680e-06,
0.334679340401e-07,
0.190904054379e-08,
0.811187736448e-10,
0.250665552389e-11,
0.546594403180e-13,
0.809426189344e-15,
0.774238295702e-17,
0.447098436530e-19,
0.141709359957e-21,
0.213765830835e-24,
0.121524412340e-27,
0.167380166790e-31,
0.183379404857e-36,
// weights 60
0.110958724796e-44,
0.243974758810e-39,
0.377162672698e-35,
0.133255961176e-31,
0.171557314767e-28,
0.102940599693e-25,
0.334575695574e-23,
0.651256725748e-21,
0.815364047300e-19,
0.692324790956e-17,
0.415244410968e-15,
0.181662457614e-13,
0.594843051597e-12,
0.148895734905e-10,
0.289935901280e-09,
0.445682277521e-08,
0.547555461926e-07,
0.543351613419e-06,
0.439428693625e-05,
0.291874190415e-04,
0.160277334681e-03,
0.731773556963e-03,
0.279132482894e-02,
0.893217836028e-02,
0.240612727660e-01,
0.547189709320e-01,
0.105298763697e+00,
0.171776156918e+00,
0.237868904958e+00,
0.279853117522e+00,
0.279853117522e+00,
0.237868904958e+00,
0.171776156918e+00,
0.105298763697e+00,
0.547189709320e-01,
0.240612727660e-01,
0.893217836028e-02,
0.279132482894e-02,
0.731773556963e-03,
0.160277334681e-03,
0.291874190415e-04,
0.439428693625e-05,
0.543351613419e-06,
0.547555461926e-07,
0.445682277521e-08,
0.289935901280e-09,
0.148895734905e-10,
0.594843051597e-12,
0.181662457614e-13,
0.415244410968e-15,
0.692324790956e-17,
0.815364047300e-19,
0.651256725748e-21,
0.334575695574e-23,
0.102940599693e-25,
0.171557314767e-28,
0.133255961176e-31,
0.377162672698e-35,
0.243974758810e-39,
0.110958724796e-44,
// weights 64
0.553570653584e-48,
0.167974799010e-42,
0.342113801099e-38,
0.155739062462e-34,
0.254966089910e-31,
0.192910359546e-28,
0.786179778889e-26,
0.191170688329e-23,
0.298286278427e-21,
0.315225456649e-19,
0.235188471067e-17,
0.128009339117e-15,
0.521862372645e-14,
0.162834073070e-12,
0.395917776693e-11,
0.761521725012e-10,
0.117361674232e-08,
0.146512531647e-07,
0.149553293672e-06,
0.125834025103e-05,
0.878849923082e-05,
0.512592913577e-04,
0.250983698512e-03,
0.103632909950e-02,
0.362258697852e-02,
0.107560405098e-01,
0.272031289536e-01,
0.587399819634e-01,
0.108498349306e+00,
0.171685842349e+00,
0.232994786062e+00,
0.271377424940e+00,
0.271377424940e+00,
0.232994786062e+00,
0.171685842349e+00,
0.108498349306e+00,
0.587399819634e-01,
0.272031289536e-01,
0.107560405098e-01,
0.362258697852e-02,
0.103632909950e-02,
0.250983698512e-03,
0.512592913577e-04,
0.878849923082e-05,
0.125834025103e-05,
0.149553293672e-06,
0.146512531647e-07,
0.117361674232e-08,
0.761521725012e-10,
0.395917776693e-11,
0.162834073070e-12,
0.521862372645e-14,
0.128009339117e-15,
0.235188471067e-17,
0.315225456649e-19,
0.298286278427e-21,
0.191170688329e-23,
0.786179778889e-26,
0.192910359546e-28,
0.254966089910e-31,
0.155739062462e-34,
0.342113801099e-38,
0.167974799010e-42,
0.553570653584e-48
};
// Points, starting with the first level
static const double gh_points[] = {
// points 1
0.0,
// points 2
-0.707106781186547524400844362105e+00,
0.707106781186547524400844362105e+00,
// points 3
-0.122474487139158904909864203735e+01,
0.0e+00,
0.122474487139158904909864203735e+01,
// points 4
-0.165068012388578455588334111112e+01,
-0.524647623275290317884060253835e+00,
0.524647623275290317884060253835e+00,
0.165068012388578455588334111112e+01,
// points 5
-0.202018287045608563292872408814e+01,
-0.958572464613818507112770593893e+00,
0.0e+00,
0.958572464613818507112770593893e+00,
0.202018287045608563292872408814e+01,
// points 6
-0.235060497367449222283392198706e+01,
-0.133584907401369694971489528297e+01,
-0.436077411927616508679215948251e+00,
0.436077411927616508679215948251e+00,
0.133584907401369694971489528297e+01,
0.235060497367449222283392198706e+01,
// points 7
-0.265196135683523349244708200652e+01,
-0.167355162876747144503180139830e+01,
-0.816287882858964663038710959027e+00,
0.0e+00,
0.816287882858964663038710959027e+00,
0.167355162876747144503180139830e+01,
0.265196135683523349244708200652e+01,
// points 8
-0.293063742025724401922350270524e+01,
-0.198165675669584292585463063977e+01,
-0.115719371244678019472076577906e+01,
-0.381186990207322116854718885584e+00,
0.381186990207322116854718885584e+00,
0.115719371244678019472076577906e+01,
0.198165675669584292585463063977e+01,
0.293063742025724401922350270524e+01,
// points 9
-0.319099320178152760723004779538e+01,
-0.226658058453184311180209693284e+01,
-0.146855328921666793166701573925e+01,
-0.723551018752837573322639864579e+00,
0.0e+00,
0.723551018752837573322639864579e+00,
0.146855328921666793166701573925e+01,
0.226658058453184311180209693284e+01,
0.319099320178152760723004779538e+01,
// points 10
-0.343615911883773760332672549432e+01,
-0.253273167423278979640896079775e+01,
-0.175668364929988177345140122011e+01,
-0.103661082978951365417749191676e+01,
-0.342901327223704608789165025557e+00,
0.342901327223704608789165025557e+00,
0.103661082978951365417749191676e+01,
0.175668364929988177345140122011e+01,
0.253273167423278979640896079775e+01,
0.343615911883773760332672549432e+01,
// points 11
-0.366847084655958251845837146485e+01,
-0.278329009978165177083671870152e+01,
-0.202594801582575533516591283121e+01,
-0.132655708449493285594973473558e+01,
-0.656809566882099765024611575383e+00,
0.0e+00,
0.656809566882099765024611575383e+00,
0.132655708449493285594973473558e+01,
0.202594801582575533516591283121e+01,
0.278329009978165177083671870152e+01,
0.366847084655958251845837146485e+01,
// points 12
-0.388972489786978191927164274724e+01,
-0.302063702512088977171067937518e+01,
-0.227950708050105990018772856942e+01,
-0.159768263515260479670966277090e+01,
-0.947788391240163743704578131060e+00,
-0.314240376254359111276611634095e+00,
0.314240376254359111276611634095e+00,
0.947788391240163743704578131060e+00,
0.159768263515260479670966277090e+01,
0.227950708050105990018772856942e+01,
0.302063702512088977171067937518e+01,
0.388972489786978191927164274724e+01,
// points 13
-0.410133759617863964117891508007e+01,
-0.324660897837240998812205115236e+01,
-0.251973568567823788343040913628e+01,
-0.185310765160151214200350644316e+01,
-0.122005503659074842622205526637e+01,
-0.605763879171060113080537108602e+00,
0.0e+00,
0.605763879171060113080537108602e+00,
0.122005503659074842622205526637e+01,
0.185310765160151214200350644316e+01,
0.251973568567823788343040913628e+01,
0.324660897837240998812205115236e+01,
0.410133759617863964117891508007e+01,
// points 14
-0.430444857047363181262129810037e+01,
-0.346265693360227055020891736115e+01,
-0.274847072498540256862499852415e+01,
-0.209518325850771681573497272630e+01,
-0.147668273114114087058350654421e+01,
-0.878713787329399416114679311861e+00,
-0.291745510672562078446113075799e+00,
0.291745510672562078446113075799e+00,
0.878713787329399416114679311861e+00,
0.147668273114114087058350654421e+01,
0.209518325850771681573497272630e+01,
0.274847072498540256862499852415e+01,
0.346265693360227055020891736115e+01,
0.430444857047363181262129810037e+01,
// points 15
-0.449999070730939155366438053053e+01,
-0.366995037340445253472922383312e+01,
-0.296716692790560324848896036355e+01,
-0.232573248617385774545404479449e+01,
-0.171999257518648893241583152515e+01,
-0.113611558521092066631913490556e+01,
-0.565069583255575748526020337198e+00,
0.0e+00,
0.565069583255575748526020337198e+00,
0.113611558521092066631913490556e+01,
0.171999257518648893241583152515e+01,
0.232573248617385774545404479449e+01,
0.296716692790560324848896036355e+01,
0.366995037340445253472922383312e+01,
0.449999070730939155366438053053e+01,
// points 16
-0.468873893930581836468849864875e+01,
-0.386944790486012269871942409801e+01,
-0.317699916197995602681399455926e+01,
-0.254620215784748136215932870545e+01,
-0.195178799091625397743465541496e+01,
-0.138025853919888079637208966969e+01,
-0.822951449144655892582454496734e+00,
-0.273481046138152452158280401965e+00,
0.273481046138152452158280401965e+00,
0.822951449144655892582454496734e+00,
0.138025853919888079637208966969e+01,
0.195178799091625397743465541496e+01,
0.254620215784748136215932870545e+01,
0.317699916197995602681399455926e+01,
0.386944790486012269871942409801e+01,
0.468873893930581836468849864875e+01,
// points 17
-0.487134519367440308834927655662e+01,
-0.406194667587547430689245559698e+01,
-0.337893209114149408338327069289e+01,
-0.275776291570388873092640349574e+01,
-0.217350282666662081927537907149e+01,
-0.161292431422123133311288254454e+01,
-0.106764872574345055363045773799e+01,
-0.531633001342654731349086553718e+00,
0.0e+00,
0.531633001342654731349086553718e+00,
0.106764872574345055363045773799e+01,
0.161292431422123133311288254454e+01,
0.217350282666662081927537907149e+01,
0.275776291570388873092640349574e+01,
0.337893209114149408338327069289e+01,
0.406194667587547430689245559698e+01,
0.487134519367440308834927655662e+01,
// points 18
-0.504836400887446676837203757885e+01,
-0.424811787356812646302342016090e+01,
-0.357376906848626607950067599377e+01,
-0.296137750553160684477863254906e+01,
-0.238629908916668600026459301424e+01,
-0.183553160426162889225383944409e+01,
-0.130092085838961736566626555439e+01,
-0.776682919267411661316659462284e+00,
-0.258267750519096759258116098711e+00,
0.258267750519096759258116098711e+00,
0.776682919267411661316659462284e+00,
0.130092085838961736566626555439e+01,
0.183553160426162889225383944409e+01,
0.238629908916668600026459301424e+01,
0.296137750553160684477863254906e+01,
0.357376906848626607950067599377e+01,
0.424811787356812646302342016090e+01,
0.504836400887446676837203757885e+01,
// points 19
-0.522027169053748216460967142500e+01,
-0.442853280660377943723498532226e+01,
-0.376218735196402009751489394104e+01,
-0.315784881834760228184318034120e+01,
-0.259113378979454256492128084112e+01,
-0.204923170985061937575050838669e+01,
-0.152417061939353303183354859367e+01,
-0.101036838713431135136859873726e+01,
-0.503520163423888209373811765050e+00,
0.0e+00,
0.503520163423888209373811765050e+00,
0.101036838713431135136859873726e+01,
0.152417061939353303183354859367e+01,
0.204923170985061937575050838669e+01,
0.259113378979454256492128084112e+01,
0.315784881834760228184318034120e+01,
0.376218735196402009751489394104e+01,
0.442853280660377943723498532226e+01,
0.522027169053748216460967142500e+01,
// points 20
-0.538748089001123286201690041068e+01,
-0.460368244955074427307767524898e+01,
-0.394476404011562521037562880052e+01,
-0.334785456738321632691492452300e+01,
-0.278880605842813048052503375640e+01,
-0.225497400208927552308233334473e+01,
-0.173853771211658620678086566214e+01,
-0.123407621539532300788581834696e+01,
-0.737473728545394358705605144252e+00,
-0.245340708300901249903836530634e+00,
0.245340708300901249903836530634e+00,
0.737473728545394358705605144252e+00,
0.123407621539532300788581834696e+01,
0.173853771211658620678086566214e+01,
0.225497400208927552308233334473e+01,
0.278880605842813048052503375640e+01,
0.334785456738321632691492452300e+01,
0.394476404011562521037562880052e+01,
0.460368244955074427307767524898e+01,
0.538748089001123286201690041068e+01,
// points 30
-6.86334529352989158106110835756e+00,
-6.13827922012393462039499237854e+00,
-5.53314715156749572511833355558e+00,
-4.98891896858994394448649710633e+00,
-4.48305535709251834188703761971e+00,
-4.00390860386122881522787601332e+00,
-3.54444387315534988692540090217e+00,
-3.09997052958644174868873332237e+00,
-2.66713212453561720057110646422e+00,
-2.24339146776150407247297999483e+00,
-1.82674114360368803883588048351e+00,
-1.41552780019818851194072510555e+00,
-1.00833827104672346180498960870e+00,
-0.603921058625552307778155678757e+00,
-0.201128576548871485545763013244e+00,
0.201128576548871485545763013244e+00,
0.603921058625552307778155678757e+00,
1.00833827104672346180498960870e+00,
1.41552780019818851194072510555e+00,
1.82674114360368803883588048351e+00,
2.24339146776150407247297999483e+00,
2.66713212453561720057110646422e+00,
3.09997052958644174868873332237e+00,
3.54444387315534988692540090217e+00,
4.00390860386122881522787601332e+00,
4.48305535709251834188703761971e+00,
4.98891896858994394448649710633e+00,
5.53314715156749572511833355558e+00,
6.13827922012393462039499237854e+00,
6.86334529352989158106110835756e+00,
// points 32
-7.12581390983e+00,
-6.40949814927e+00,
-5.81222594952e+00,
-5.27555098652e+00,
-4.77716450350e+00,
-4.30554795335e+00,
-3.85375548547e+00,
-3.41716749282e+00,
-2.99249082500e+00,
-2.57724953773e+00,
-2.16949918361e+00,
-1.76765410946e+00,
-1.37037641095e+00,
-0.976500463590e+00,
-0.584978765436e+00,
-0.194840741569e+00,
0.194840741569e+00,
0.584978765436e+00,
0.976500463590e+00,
1.37037641095e+00,
1.76765410946e+00,
2.16949918361e+00,
2.57724953773e+00,
2.99249082500e+00,
3.41716749282e+00,
3.85375548547e+00,
4.30554795335e+00,
4.77716450350e+00,
5.27555098652e+00,
5.81222594952e+00,
6.40949814927e+00,
7.12581390983e+00,
// points 40
-8.09876113925e+00,
-7.41158253149e+00,
-6.84023730525e+00,
-6.32825535122e+00,
-5.85409505603e+00,
-5.40665424797e+00,
-4.97926097855e+00,
-4.56750207284e+00,
-4.16825706683e+00,
-3.77920675344e+00,
-3.39855826586e+00,
-3.02487988390e+00,
-2.65699599844e+00,
-2.29391714188e+00,
-1.93479147228e+00,
-1.57886989493e+00,
-1.22548010905e+00,
-0.874006612357e+00,
-0.523874713832e+00,
-0.174537214598e+00,
0.174537214598e+00,
0.523874713832e+00,
0.874006612357e+00,
1.22548010905e+00,
1.57886989493e+00,
1.93479147228e+00,
2.29391714188e+00,
2.65699599844e+00,
3.02487988390e+00,
3.39855826586e+00,
3.77920675344e+00,
4.16825706683e+00,
4.56750207284e+00,
4.97926097855e+00,
5.40665424797e+00,
5.85409505603e+00,
6.32825535122e+00,
6.84023730525e+00,
7.41158253149e+00,
8.09876113925e+00,
// points 50
-9.18240695813e+00,
-8.52277103092e+00,
-7.97562236821e+00,
-7.48640942986e+00,
-7.03432350977e+00,
-6.60864797386e+00,
-6.20295251927e+00,
-5.81299467542e+00,
-5.43578608722e+00,
-5.06911758492e+00,
-4.71129366617e+00,
-4.36097316045e+00,
-4.01706817286e+00,
-3.67867706252e+00,
-3.34503831394e+00,
-3.01549776957e+00,
-2.68948470227e+00,
-2.36649390430e+00,
-2.04607196869e+00,
-1.72780654752e+00,
-1.41131775490e+00,
-1.09625112896e+00,
-0.782271729555e+00,
-0.469059056678e+00,
-0.156302546889e+00,
0.156302546889e+00,
0.469059056678e+00,
0.782271729555e+00,
1.09625112896e+00,
1.41131775490e+00,
1.72780654752e+00,
2.04607196869e+00,
2.36649390430e+00,
2.68948470227e+00,
3.01549776957e+00,
3.34503831394e+00,
3.67867706252e+00,
4.01706817286e+00,
4.36097316045e+00,
4.71129366617e+00,
5.06911758492e+00,
5.43578608722e+00,
5.81299467542e+00,
6.20295251927e+00,
6.60864797386e+00,
7.03432350977e+00,
7.48640942986e+00,
7.97562236821e+00,
8.52277103092e+00,
9.18240695813e+00,
// points 60
-10.1591092462e+00,
-9.52090367701e+00,
-8.99239800140e+00,
-8.52056928412e+00,
-8.08518865425e+00,
-7.67583993750e+00,
-7.28627659440e+00,
-6.91238153219e+00,
-6.55125916706e+00,
-6.20077355799e+00,
-5.85929019639e+00,
-5.52552108614e+00,
-5.19842653458e+00,
-4.87715007747e+00,
-4.56097375794e+00,
-4.24928643596e+00,
-3.94156073393e+00,
-3.63733587617e+00,
-3.33620465355e+00,
-3.03780333823e+00,
-2.74180374807e+00,
-2.44790690231e+00,
-2.15583787123e+00,
-1.86534153123e+00,
-1.57617901198e+00,
-1.28812467487e+00,
-1.00096349956e+00,
-0.714488781673e+00,
-0.428500064221e+00,
-0.142801238703e+00,
0.142801238703e+00,
0.428500064221e+00,
0.714488781673e+00,
1.00096349956e+00,
1.28812467487e+00,
1.57617901198e+00,
1.86534153123e+00,
2.15583787123e+00,
2.44790690231e+00,
2.74180374807e+00,
3.03780333823e+00,
3.33620465355e+00,
3.63733587617e+00,
3.94156073393e+00,
4.24928643596e+00,
4.56097375794e+00,
4.87715007747e+00,
5.19842653458e+00,
5.52552108614e+00,
5.85929019639e+00,
6.20077355799e+00,
6.55125916706e+00,
6.91238153219e+00,
7.28627659440e+00,
7.67583993750e+00,
8.08518865425e+00,
8.52056928412e+00,
8.99239800140e+00,
9.52090367701e+00,
10.1591092462e+00,
// points 64
-10.5261231680e+00,
-9.89528758683e+00,
-9.37315954965e+00,
-8.90724909996e+00,
-8.47752908338e+00,
-8.07368728501e+00,
-7.68954016404e+00,
-7.32101303278e+00,
-6.96524112055e+00,
-6.62011226264e+00,
-6.28401122877e+00,
-5.95566632680e+00,
-5.63405216435e+00,
-5.31832522463e+00,
-5.00777960220e+00,
-4.70181564741e+00,
-4.39991716823e+00,
-4.10163447457e+00,
-3.80657151395e+00,
-3.51437593574e+00,
-3.22473129199e+00,
-2.93735082300e+00,
-2.65197243543e+00,
-2.36835458863e+00,
-2.08627287988e+00,
-1.80551717147e+00,
-1.52588914021e+00,
-1.24720015694e+00,
-0.969269423071e+00,
-0.691922305810e+00,
-0.414988824121e+00,
-0.138302244987e+00,
0.138302244987e+00,
0.414988824121e+00,
0.691922305810e+00,
0.969269423071e+00,
1.24720015694e+00,
1.52588914021e+00,
1.80551717147e+00,
2.08627287988e+00,
2.36835458863e+00,
2.65197243543e+00,
2.93735082300e+00,
3.22473129199e+00,
3.51437593574e+00,
3.80657151395e+00,
4.10163447457e+00,
4.39991716823e+00,
4.70181564741e+00,
5.00777960220e+00,
5.31832522463e+00,
5.63405216435e+00,
5.95566632680e+00,
6.28401122877e+00,
6.62011226264e+00,
6.96524112055e+00,
7.32101303278e+00,
7.68954016404e+00,
8.07368728501e+00,
8.47752908338e+00,
8.90724909996e+00,
9.37315954965e+00,
9.89528758683e+00,
10.5261231680e+00
};
// Gauss-Legendre quadrature; prefix gl
// Number of levels
static const int gl_num_levels = 22;
// Number of points in each level
static const int gl_num_points[] = {
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
32, 64
};
// Weights, starting with the first level
static const double gl_weights[] = {
// weight 1
2.0e+00,
// weights 2
1.0e+00,
1.0e+00,
// weights 3
0.555555555555555555555555555555e+00,
0.888888888888888888888888888888e+00,
0.555555555555555555555555555555e+00,
// weights 4
0.347854845137453857373063949222e+00,
0.652145154862546142626936050778e+00,
0.652145154862546142626936050778e+00,
0.347854845137453857373063949222e+00,
// weights 5
0.236926885056189087514264040720e+00,
0.478628670499366468041291514836e+00,
0.568888888888888888888888888889e+00,
0.478628670499366468041291514836e+00,
0.236926885056189087514264040720e+00,
// weights 6
0.171324492379170345040296142173e+00,
0.360761573048138607569833513838e+00,
0.467913934572691047389870343990e+00,
0.467913934572691047389870343990e+00,
0.360761573048138607569833513838e+00,
0.171324492379170345040296142173e+00,
// weights 7
0.129484966168869693270611432679e+00,
0.279705391489276667901467771424e+00,
0.381830050505118944950369775489e+00,
0.417959183673469387755102040816e+00,
0.381830050505118944950369775489e+00,
0.279705391489276667901467771424e+00,
0.129484966168869693270611432679e+00,
// weights 8
0.101228536290376259152531354310e+00,
0.222381034453374470544355994426e+00,
0.313706645877887287337962201987e+00,
0.362683783378361982965150449277e+00,
0.362683783378361982965150449277e+00,
0.313706645877887287337962201987e+00,
0.222381034453374470544355994426e+00,
0.101228536290376259152531354310e+00,
// weights 9
0.812743883615744119718921581105e-01,
0.180648160694857404058472031243e+00,
0.260610696402935462318742869419e+00,
0.312347077040002840068630406584e+00,
0.330239355001259763164525069287e+00,
0.312347077040002840068630406584e+00,
0.260610696402935462318742869419e+00,
0.180648160694857404058472031243e+00,
0.812743883615744119718921581105e-01,
// weights 10
0.666713443086881375935688098933e-01,
0.149451349150580593145776339658e+00,
0.219086362515982043995534934228e+00,
0.269266719309996355091226921569e+00,
0.295524224714752870173892994651e+00,
0.295524224714752870173892994651e+00,
0.269266719309996355091226921569e+00,
0.219086362515982043995534934228e+00,
0.149451349150580593145776339658e+00,
0.666713443086881375935688098933e-01,
// weights 11
0.556685671161736664827537204425e-01,
0.125580369464904624634694299224e+00,
0.186290210927734251426097641432e+00,
0.233193764591990479918523704843e+00,
0.262804544510246662180688869891e+00,
0.272925086777900630714483528336e+00,
0.262804544510246662180688869891e+00,
0.233193764591990479918523704843e+00,
0.186290210927734251426097641432e+00,
0.125580369464904624634694299224e+00,
0.556685671161736664827537204425e-01,
// weights 12
0.471753363865118271946159614850e-01,
0.106939325995318430960254718194e+00,
0.160078328543346226334652529543e+00,
0.203167426723065921749064455810e+00,
0.233492536538354808760849898925e+00,
0.249147045813402785000562436043e+00,
0.249147045813402785000562436043e+00,
0.233492536538354808760849898925e+00,
0.203167426723065921749064455810e+00,
0.160078328543346226334652529543e+00,
0.106939325995318430960254718194e+00,
0.471753363865118271946159614850e-01,
// weights 13
0.404840047653158795200215922010e-01,
0.921214998377284479144217759538e-01,
0.138873510219787238463601776869e+00,
0.178145980761945738280046691996e+00,
0.207816047536888502312523219306e+00,
0.226283180262897238412090186040e+00,
0.232551553230873910194589515269e+00,
0.226283180262897238412090186040e+00,
0.207816047536888502312523219306e+00,
0.178145980761945738280046691996e+00,
0.138873510219787238463601776869e+00,
0.921214998377284479144217759538e-01,
0.404840047653158795200215922010e-01,
// weights 14
0.351194603317518630318328761382e-01,
0.801580871597602098056332770629e-01,
0.121518570687903184689414809072e+00,
0.157203167158193534569601938624e+00,
0.185538397477937813741716590125e+00,
0.205198463721295603965924065661e+00,
0.215263853463157790195876443316e+00,
0.215263853463157790195876443316e+00,
0.205198463721295603965924065661e+00,
0.185538397477937813741716590125e+00,
0.157203167158193534569601938624e+00,
0.121518570687903184689414809072e+00,
0.801580871597602098056332770629e-01,
0.351194603317518630318328761382e-01,
// weights 15
0.307532419961172683546283935772e-01,
0.703660474881081247092674164507e-01,
0.107159220467171935011869546686e+00,
0.139570677926154314447804794511e+00,
0.166269205816993933553200860481e+00,
0.186161000015562211026800561866e+00,
0.198431485327111576456118326444e+00,
0.202578241925561272880620199968e+00,
0.198431485327111576456118326444e+00,
0.186161000015562211026800561866e+00,
0.166269205816993933553200860481e+00,
0.139570677926154314447804794511e+00,
0.107159220467171935011869546686e+00,
0.703660474881081247092674164507e-01,
0.307532419961172683546283935772e-01,
// weights 16
0.271524594117540948517805724560e-01,
0.622535239386478928628438369944e-01,
0.951585116824927848099251076022e-01,
0.124628971255533872052476282192e+00,
0.149595988816576732081501730547e+00,
0.169156519395002538189312079030e+00,
0.182603415044923588866763667969e+00,
0.189450610455068496285396723208e+00,
0.189450610455068496285396723208e+00,
0.182603415044923588866763667969e+00,
0.169156519395002538189312079030e+00,
0.149595988816576732081501730547e+00,
0.124628971255533872052476282192e+00,
0.951585116824927848099251076022e-01,
0.622535239386478928628438369944e-01,
0.271524594117540948517805724560e-01,
// weights 17
0.241483028685479319601100262876e-01,
0.554595293739872011294401653582e-01,
0.850361483171791808835353701911e-01,
0.111883847193403971094788385626e+00,
0.135136368468525473286319981702e+00,
0.154045761076810288081431594802e+00,
0.168004102156450044509970663788e+00,
0.176562705366992646325270990113e+00,
0.179446470356206525458265644262e+00,
0.176562705366992646325270990113e+00,
0.168004102156450044509970663788e+00,
0.154045761076810288081431594802e+00,
0.135136368468525473286319981702e+00,
0.111883847193403971094788385626e+00,
0.850361483171791808835353701911e-01,
0.554595293739872011294401653582e-01,
0.241483028685479319601100262876e-01,
// weights 18
0.216160135264833103133427102665e-01,
0.497145488949697964533349462026e-01,
0.764257302548890565291296776166e-01,
0.100942044106287165562813984925e+00,
0.122555206711478460184519126800e+00,
0.140642914670650651204731303752e+00,
0.154684675126265244925418003836e+00,
0.164276483745832722986053776466e+00,
0.169142382963143591840656470135e+00,
0.169142382963143591840656470135e+00,
0.164276483745832722986053776466e+00,
0.154684675126265244925418003836e+00,
0.140642914670650651204731303752e+00,
0.122555206711478460184519126800e+00,
0.100942044106287165562813984925e+00,
0.764257302548890565291296776166e-01,
0.497145488949697964533349462026e-01,
0.216160135264833103133427102665e-01,
// weights 19
0.194617882297264770363120414644e-01,
0.448142267656996003328381574020e-01,
0.690445427376412265807082580060e-01,
0.914900216224499994644620941238e-01,
0.111566645547333994716023901682e+00,
0.128753962539336227675515784857e+00,
0.142606702173606611775746109442e+00,
0.152766042065859666778855400898e+00,
0.158968843393954347649956439465e+00,
0.161054449848783695979163625321e+00,
0.158968843393954347649956439465e+00,
0.152766042065859666778855400898e+00,
0.142606702173606611775746109442e+00,
0.128753962539336227675515784857e+00,
0.111566645547333994716023901682e+00,
0.914900216224499994644620941238e-01,
0.690445427376412265807082580060e-01,
0.448142267656996003328381574020e-01,
0.194617882297264770363120414644e-01,
// weights 20
0.176140071391521183118619623519e-01,
0.406014298003869413310399522749e-01,
0.626720483341090635695065351870e-01,
0.832767415767047487247581432220e-01,
0.101930119817240435036750135480e+00,
0.118194531961518417312377377711e+00,
0.131688638449176626898494499748e+00,
0.142096109318382051329298325067e+00,
0.149172986472603746787828737002e+00,
0.152753387130725850698084331955e+00,
0.152753387130725850698084331955e+00,
0.149172986472603746787828737002e+00,
0.142096109318382051329298325067e+00,
0.131688638449176626898494499748e+00,
0.118194531961518417312377377711e+00,
0.101930119817240435036750135480e+00,
0.832767415767047487247581432220e-01,
0.626720483341090635695065351870e-01,
0.406014298003869413310399522749e-01,
0.176140071391521183118619623519e-01,
// weights 32
0.701861000947009660040706373885e-02,
0.162743947309056706051705622064e-01,
0.253920653092620594557525897892e-01,
0.342738629130214331026877322524e-01,
0.428358980222266806568786466061e-01,
0.509980592623761761961632446895e-01,
0.586840934785355471452836373002e-01,
0.658222227763618468376500637069e-01,
0.723457941088485062253993564785e-01,
0.781938957870703064717409188283e-01,
0.833119242269467552221990746043e-01,
0.876520930044038111427714627518e-01,
0.911738786957638847128685771116e-01,
0.938443990808045656391802376681e-01,
0.956387200792748594190820022041e-01,
0.965400885147278005667648300636e-01,
0.965400885147278005667648300636e-01,
0.956387200792748594190820022041e-01,
0.938443990808045656391802376681e-01,
0.911738786957638847128685771116e-01,
0.876520930044038111427714627518e-01,
0.833119242269467552221990746043e-01,
0.781938957870703064717409188283e-01,
0.723457941088485062253993564785e-01,
0.658222227763618468376500637069e-01,
0.586840934785355471452836373002e-01,
0.509980592623761761961632446895e-01,
0.428358980222266806568786466061e-01,
0.342738629130214331026877322524e-01,
0.253920653092620594557525897892e-01,
0.162743947309056706051705622064e-01,
0.701861000947009660040706373885e-02,
// weights 64
0.178328072169643294729607914497e-02,
0.414703326056246763528753572855e-02,
0.650445796897836285611736039998e-02,
0.884675982636394772303091465973e-02,
0.111681394601311288185904930192e-01,
0.134630478967186425980607666860e-01,
0.157260304760247193219659952975e-01,
0.179517157756973430850453020011e-01,
0.201348231535302093723403167285e-01,
0.222701738083832541592983303842e-01,
0.243527025687108733381775504091e-01,
0.263774697150546586716917926252e-01,
0.283396726142594832275113052002e-01,
0.302346570724024788679740598195e-01,
0.320579283548515535854675043479e-01,
0.338051618371416093915654821107e-01,
0.354722132568823838106931467152e-01,
0.370551285402400460404151018096e-01,
0.385501531786156291289624969468e-01,
0.399537411327203413866569261283e-01,
0.412625632426235286101562974736e-01,
0.424735151236535890073397679088e-01,
0.435837245293234533768278609737e-01,
0.445905581637565630601347100309e-01,
0.454916279274181444797709969713e-01,
0.462847965813144172959532492323e-01,
0.469681828162100173253262857546e-01,
0.475401657148303086622822069442e-01,
0.479993885964583077281261798713e-01,
0.483447622348029571697695271580e-01,
0.485754674415034269347990667840e-01,
0.486909570091397203833653907347e-01,
0.486909570091397203833653907347e-01,
0.485754674415034269347990667840e-01,
0.483447622348029571697695271580e-01,
0.479993885964583077281261798713e-01,
0.475401657148303086622822069442e-01,
0.469681828162100173253262857546e-01,
0.462847965813144172959532492323e-01,
0.454916279274181444797709969713e-01,
0.445905581637565630601347100309e-01,
0.435837245293234533768278609737e-01,
0.424735151236535890073397679088e-01,
0.412625632426235286101562974736e-01,
0.399537411327203413866569261283e-01,
0.385501531786156291289624969468e-01,
0.370551285402400460404151018096e-01,
0.354722132568823838106931467152e-01,
0.338051618371416093915654821107e-01,
0.320579283548515535854675043479e-01,
0.302346570724024788679740598195e-01,
0.283396726142594832275113052002e-01,
0.263774697150546586716917926252e-01,
0.243527025687108733381775504091e-01,
0.222701738083832541592983303842e-01,
0.201348231535302093723403167285e-01,
0.179517157756973430850453020011e-01,
0.157260304760247193219659952975e-01,
0.134630478967186425980607666860e-01,
0.111681394601311288185904930192e-01,
0.884675982636394772303091465973e-02,
0.650445796897836285611736039998e-02,
0.414703326056246763528753572855e-02,
0.178328072169643294729607914497e-02
};
// Points, starting with the first level
static const double gl_points[] = {
// points 1
0.0e+00,
// points 2
-0.577350269189625764509148780502e+00,
0.577350269189625764509148780502e+00,
// points 3
-0.774596669241483377035853079956e+00,
0.0e+00,
0.774596669241483377035853079956e+00,
// points 4
-0.861136311594052575223946488893e+00,
-0.339981043584856264802665759103e+00,
0.339981043584856264802665759103e+00,
0.861136311594052575223946488893e+00,
// points 5
-0.906179845938663992797626878299e+00,
-0.538469310105683091036314420700e+00,
0.0e+00,
0.538469310105683091036314420700e+00,
0.906179845938663992797626878299e+00,
// points 6
-0.932469514203152027812301554494e+00,
-0.661209386466264513661399595020e+00,
-0.238619186083196908630501721681e+00,
0.238619186083196908630501721681e+00,
0.661209386466264513661399595020e+00,
0.932469514203152027812301554494e+00,
// points 7
-0.949107912342758524526189684048e+00,
-0.741531185599394439863864773281e+00,
-0.405845151377397166906606412077e+00,
0.0e+00,
0.405845151377397166906606412077e+00,
0.741531185599394439863864773281e+00,
0.949107912342758524526189684048e+00,
// points 8
-0.960289856497536231683560868569e+00,
-0.796666477413626739591553936476e+00,
-0.525532409916328985817739049189e+00,
-0.183434642495649804939476142360e+00,
0.183434642495649804939476142360e+00,
0.525532409916328985817739049189e+00,
0.796666477413626739591553936476e+00,
0.960289856497536231683560868569e+00,
// points 9
-0.968160239507626089835576202904e+00,
-0.836031107326635794299429788070e+00,
-0.613371432700590397308702039341e+00,
-0.324253423403808929038538014643e+00,
0.0e+00,
0.324253423403808929038538014643e+00,
0.613371432700590397308702039341e+00,
0.836031107326635794299429788070e+00,
0.968160239507626089835576202904e+00,
// points 10
-0.973906528517171720077964012084e+00,
-0.865063366688984510732096688423e+00,
-0.679409568299024406234327365115e+00,
-0.433395394129247190799265943166e+00,
-0.148874338981631210884826001130e+00,
0.148874338981631210884826001130e+00,
0.433395394129247190799265943166e+00,
0.679409568299024406234327365115e+00,
0.865063366688984510732096688423e+00,
0.973906528517171720077964012084e+00,
// points 11
-0.978228658146056992803938001123e+00,
-0.887062599768095299075157769304e+00,
-0.730152005574049324093416252031e+00,
-0.519096129206811815925725669459e+00,
-0.269543155952344972331531985401e+00,
0.0e+00,
0.269543155952344972331531985401e+00,
0.519096129206811815925725669459e+00,
0.730152005574049324093416252031e+00,
0.887062599768095299075157769304e+00,
0.978228658146056992803938001123e+00,
// points 12
-0.981560634246719250690549090149e+00,
-0.904117256370474856678465866119e+00,
-0.769902674194304687036893833213e+00,
-0.587317954286617447296702418941e+00,
-0.367831498998180193752691536644e+00,
-0.125233408511468915472441369464e+00,
0.125233408511468915472441369464e+00,
0.367831498998180193752691536644e+00,
0.587317954286617447296702418941e+00,
0.769902674194304687036893833213e+00,
0.904117256370474856678465866119e+00,
0.981560634246719250690549090149e+00,
// points 13
-0.984183054718588149472829448807e+00,
-0.917598399222977965206547836501e+00,
-0.801578090733309912794206489583e+00,
-0.642349339440340220643984606996e+00,
-0.448492751036446852877912852128e+00,
-0.230458315955134794065528121098e+00,
0.0e+00,
0.230458315955134794065528121098e+00,
0.448492751036446852877912852128e+00,
0.642349339440340220643984606996e+00,
0.801578090733309912794206489583e+00,
0.917598399222977965206547836501e+00,
0.984183054718588149472829448807e+00,
// points 14
-0.986283808696812338841597266704e+00,
-0.928434883663573517336391139378e+00,
-0.827201315069764993189794742650e+00,
-0.687292904811685470148019803019e+00,
-0.515248636358154091965290718551e+00,
-0.319112368927889760435671824168e+00,
-0.108054948707343662066244650220e+00,
0.108054948707343662066244650220e+00,
0.319112368927889760435671824168e+00,
0.515248636358154091965290718551e+00,
0.687292904811685470148019803019e+00,
0.827201315069764993189794742650e+00,
0.928434883663573517336391139378e+00,
0.986283808696812338841597266704e+00,
// points 15
-0.987992518020485428489565718587e+00,
-0.937273392400705904307758947710e+00,
-0.848206583410427216200648320774e+00,
-0.724417731360170047416186054614e+00,
-0.570972172608538847537226737254e+00,
-0.394151347077563369897207370981e+00,
-0.201194093997434522300628303395e+00,
0.0e+00,
0.201194093997434522300628303395e+00,
0.394151347077563369897207370981e+00,
0.570972172608538847537226737254e+00,
0.724417731360170047416186054614e+00,
0.848206583410427216200648320774e+00,
0.937273392400705904307758947710e+00,
0.987992518020485428489565718587e+00,
// points 16
-0.989400934991649932596154173450e+00,
-0.944575023073232576077988415535e+00,
-0.865631202387831743880467897712e+00,
-0.755404408355003033895101194847e+00,
-0.617876244402643748446671764049e+00,
-0.458016777657227386342419442984e+00,
-0.281603550779258913230460501460e+00,
-0.950125098376374401853193354250e-01,
0.950125098376374401853193354250e-01,
0.281603550779258913230460501460e+00,
0.458016777657227386342419442984e+00,
0.617876244402643748446671764049e+00,
0.755404408355003033895101194847e+00,
0.865631202387831743880467897712e+00,
0.944575023073232576077988415535e+00,
0.989400934991649932596154173450e+00,
// points 17
-0.990575475314417335675434019941e+00,
-0.950675521768767761222716957896e+00,
-0.880239153726985902122955694488e+00,
-0.781514003896801406925230055520e+00,
-0.657671159216690765850302216643e+00,
-0.512690537086476967886246568630e+00,
-0.351231763453876315297185517095e+00,
-0.178484181495847855850677493654e+00,
0.0e+00,
0.178484181495847855850677493654e+00,
0.351231763453876315297185517095e+00,
0.512690537086476967886246568630e+00,
0.657671159216690765850302216643e+00,
0.781514003896801406925230055520e+00,
0.880239153726985902122955694488e+00,
0.950675521768767761222716957896e+00,
0.990575475314417335675434019941e+00,
// points 18
-0.991565168420930946730016004706e+00,
-0.955823949571397755181195892930e+00,
-0.892602466497555739206060591127e+00,
-0.803704958972523115682417455015e+00,
-0.691687043060353207874891081289e+00,
-0.559770831073947534607871548525e+00,
-0.411751161462842646035931793833e+00,
-0.251886225691505509588972854878e+00,
-0.847750130417353012422618529358e-01,
0.847750130417353012422618529358e-01,
0.251886225691505509588972854878e+00,
0.411751161462842646035931793833e+00,
0.559770831073947534607871548525e+00,
0.691687043060353207874891081289e+00,
0.803704958972523115682417455015e+00,
0.892602466497555739206060591127e+00,
0.955823949571397755181195892930e+00,
0.991565168420930946730016004706e+00,
// points 19
-0.992406843843584403189017670253e+00,
-0.960208152134830030852778840688e+00,
-0.903155903614817901642660928532e+00,
-0.822714656537142824978922486713e+00,
-0.720966177335229378617095860824e+00,
-0.600545304661681023469638164946e+00,
-0.464570741375960945717267148104e+00,
-0.316564099963629831990117328850e+00,
-0.160358645640225375868096115741e+00,
0.0e+00,
0.160358645640225375868096115741e+00,
0.316564099963629831990117328850e+00,
0.464570741375960945717267148104e+00,
0.600545304661681023469638164946e+00,
0.720966177335229378617095860824e+00,
0.822714656537142824978922486713e+00,
0.903155903614817901642660928532e+00,
0.960208152134830030852778840688e+00,
0.992406843843584403189017670253e+00,
// points 20
-0.993128599185094924786122388471e+00,
-0.963971927277913791267666131197e+00,
-0.912234428251325905867752441203e+00,
-0.839116971822218823394529061702e+00,
-0.746331906460150792614305070356e+00,
-0.636053680726515025452836696226e+00,
-0.510867001950827098004364050955e+00,
-0.373706088715419560672548177025e+00,
-0.227785851141645078080496195369e+00,
-0.765265211334973337546404093988e-01,
0.765265211334973337546404093988e-01,
0.227785851141645078080496195369e+00,
0.373706088715419560672548177025e+00,
0.510867001950827098004364050955e+00,
0.636053680726515025452836696226e+00,
0.746331906460150792614305070356e+00,
0.839116971822218823394529061702e+00,
0.912234428251325905867752441203e+00,
0.963971927277913791267666131197e+00,
0.993128599185094924786122388471e+00,
// points 32
-0.997263861849481563544981128665e+00,
-0.985611511545268335400175044631e+00,
-0.964762255587506430773811928118e+00,
-0.934906075937739689170919134835e+00,
-0.896321155766052123965307243719e+00,
-0.849367613732569970133693004968e+00,
-0.794483795967942406963097298970e+00,
-0.732182118740289680387426665091e+00,
-0.663044266930215200975115168663e+00,
-0.587715757240762329040745476402e+00,
-0.506899908932229390023747474378e+00,
-0.421351276130635345364119436172e+00,
-0.331868602282127649779916805730e+00,
-0.239287362252137074544603209166e+00,
-0.144471961582796493485186373599e+00,
-0.483076656877383162348125704405e-01,
0.483076656877383162348125704405e-01,
0.144471961582796493485186373599e+00,
0.239287362252137074544603209166e+00,
0.331868602282127649779916805730e+00,
0.421351276130635345364119436172e+00,
0.506899908932229390023747474378e+00,
0.587715757240762329040745476402e+00,
0.663044266930215200975115168663e+00,
0.732182118740289680387426665091e+00,
0.794483795967942406963097298970e+00,
0.849367613732569970133693004968e+00,
0.896321155766052123965307243719e+00,
0.934906075937739689170919134835e+00,
0.964762255587506430773811928118e+00,
0.985611511545268335400175044631e+00,
0.997263861849481563544981128665e+00,
// points 64
-0.999305041735772139456905624346e+00,
-0.996340116771955279346924500676e+00,
-0.991013371476744320739382383443e+00,
-0.983336253884625956931299302157e+00,
-0.973326827789910963741853507352e+00,
-0.961008799652053718918614121897e+00,
-0.946411374858402816062481491347e+00,
-0.929569172131939575821490154559e+00,
-0.910522137078502805756380668008e+00,
-0.889315445995114105853404038273e+00,
-0.865999398154092819760783385070e+00,
-0.840629296252580362751691544696e+00,
-0.813265315122797559741923338086e+00,
-0.783972358943341407610220525214e+00,
-0.752819907260531896611863774886e+00,
-0.719881850171610826848940217832e+00,
-0.685236313054233242563558371031e+00,
-0.648965471254657339857761231993e+00,
-0.611155355172393250248852971019e+00,
-0.571895646202634034283878116659e+00,
-0.531279464019894545658013903544e+00,
-0.489403145707052957478526307022e+00,
-0.446366017253464087984947714759e+00,
-0.402270157963991603695766771260e+00,
-0.357220158337668115950442615046e+00,
-0.311322871990210956157512698560e+00,
-0.264687162208767416373964172510e+00,
-0.217423643740007084149648748989e+00,
-0.169644420423992818037313629748e+00,
-0.121462819296120554470376463492e+00,
-0.729931217877990394495429419403e-01,
-0.243502926634244325089558428537e-01,
0.243502926634244325089558428537e-01,
0.729931217877990394495429419403e-01,
0.121462819296120554470376463492e+00,
0.169644420423992818037313629748e+00,
0.217423643740007084149648748989e+00,
0.264687162208767416373964172510e+00,
0.311322871990210956157512698560e+00,
0.357220158337668115950442615046e+00,
0.402270157963991603695766771260e+00,
0.446366017253464087984947714759e+00,
0.489403145707052957478526307022e+00,
0.531279464019894545658013903544e+00,
0.571895646202634034283878116659e+00,
0.611155355172393250248852971019e+00,
0.648965471254657339857761231993e+00,
0.685236313054233242563558371031e+00,
0.719881850171610826848940217832e+00,
0.752819907260531896611863774886e+00,
0.783972358943341407610220525214e+00,
0.813265315122797559741923338086e+00,
0.840629296252580362751691544696e+00,
0.865999398154092819760783385070e+00,
0.889315445995114105853404038273e+00,
0.910522137078502805756380668008e+00,
0.929569172131939575821490154559e+00,
0.946411374858402816062481491347e+00,
0.961008799652053718918614121897e+00,
0.973326827789910963741853507352e+00,
0.983336253884625956931299302157e+00,
0.991013371476744320739382383443e+00,
0.996340116771955279346924500676e+00,
0.999305041735772139456905624346e+00
};
dynare++/integ/cc/product.cc deleted 100644 → 0
View file @ 53ae5493
/*
* Copyright © 2005 Ondra Kamenik
* Copyright © 2019 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/>.
*/
#include "product.hh"
#include "symmetry.hh"
#include <iostream>
#include <iomanip>
/* This constructs a product iterator corresponding to index (j0,0,…,0). */
prodpit::prodpit(const ProductQuadrature &q, int j0, int l)
: prodq(q), level(l), npoints(q.uquad.numPoints(l)),
jseq(q.dimen(), 0),
end_flag(false),
sig{q.dimen()},
p{q.dimen()}
{
if (j0 < npoints)
{
jseq[0] = j0;
setPointAndWeight();
}
else
end_flag = true;
}
bool
prodpit::operator==(const prodpit &ppit) const
{
return &prodq == &ppit.prodq && end_flag == ppit.end_flag && jseq == ppit.jseq;
}
prodpit &
prodpit::operator++()
{
int i = prodq.dimen()-1;
jseq[i]++;
while (i >= 0 && jseq[i] == npoints)
{
jseq[i] = 0;
i--;
if (i >= 0)
jseq[i]++;
}
sig.signalAfter(std::max(i, 0));
if (i == -1)
end_flag = true;
if (!end_flag)
setPointAndWeight();
return *this;
}
/* This calculates the weight and sets point coordinates from the indices. */
void
prodpit::setPointAndWeight()
{
w = 1.0;
for (int i = 0; i < prodq.dimen(); i++)
{
p[i] = (prodq.uquad).point(level, jseq[i]);
w *= (prodq.uquad).weight(level, jseq[i]);
}
}
/* Debug print. */
void
prodpit::print() const
{
auto ff = std::cout.flags();
std::cout << "j=[";
for (int i = 0; i < prodq.dimen(); i++)
std::cout << std::setw(2) << jseq[i];
std::cout << std::showpos << std::fixed << std::setprecision(3)
<< "] " << std::setw(4) << w << "*(";
for (int i = 0; i < prodq.dimen()-1; i++)
std::cout << std::setw(4) << p[i] << ' ';
std::cout << std::setw(4) << p[prodq.dimen()-1] << ')' << std::endl;
std::cout.flags(ff);
}
ProductQuadrature::ProductQuadrature(int d, const OneDQuadrature &uq)
: QuadratureImpl<prodpit>(d), uquad(uq)
{
// TODO: check d≥1
}
/* This calls prodpit constructor to return an iterator which points
approximatelly at ‘ti’-th portion out of ‘tn’ portions. First we find
out how many points are in the level, and then construct an interator
(j0,0,…,0) where j0=ti·npoints/tn. */
prodpit
ProductQuadrature::begin(int ti, int tn, int l) const
{
// TODO: raise if l<dimen()
// TODO: check l ≤ uquad.numLevels()
int npoints = uquad.numPoints(l);
return prodpit(*this, ti*npoints/tn, l);
}
/* This just starts at the first level and goes to a higher level as long as a
number of evaluations (which is nₖᵈ for k being the level) is less than the
given number of evaluations. */
void
ProductQuadrature::designLevelForEvals(int max_evals, int &lev, int &evals) const
{
int last_evals;
evals = 1;
lev = 1;
do
{
lev++;
last_evals = evals;
evals = numEvals(lev);
}
while (lev < uquad.numLevels()-2 && evals < max_evals);
lev--;
evals = last_evals;
}
dynare++/integ/cc/product.hh deleted 100644 → 0
View file @ 53ae5493
/*
* Copyright © 2005 Ondra Kamenik
* Copyright © 2019 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/>.
*/
// Product quadrature.
/* This file defines a product multidimensional quadrature. If Qₖ$ denotes the
one dimensional quadrature, then the product quadrature Q of k level and
dimension d takes the form:
nₖ nₖ
Qf = ∑ … ∑ w_i₁·…·w_{i_d} f(x_i₁,…,x_{i_d})
i₁=1 i_d=1
which can be written in terms of the one dimensional quadrature Qₖ as
Qf=(Qₖ⊗…⊗Qₖ)f
Here we define the product quadrature iterator prodpit and plug it into
QuadratureImpl to obtains ProductQuadrature. */
#ifndef PRODUCT_H
#define PRODUCT_H
#include "int_sequence.hh"
#include "vector_function.hh"
#include "quadrature.hh"
/* This defines a product point iterator. We have to maintain the following: a
pointer to product quadrature in order to know the dimension and the
underlying one dimensional quadrature, then level, number of points in the
level, integer sequence of indices, signal, the coordinates of the point and
the weight.
The point indices, signal, and point coordinates are implmented as pointers
in order to allow for empty constructor.
The constructor prodpit(const ProductQuadrature& q, int j0, int l)
constructs an iterator pointing to (j0,0,…,0), which is used by begin()
dictated by QuadratureImpl. */
class ProductQuadrature;
class prodpit
{
protected:
const ProductQuadrature &prodq;
int level{0};
int npoints{0};
IntSequence jseq;
bool end_flag{true};
ParameterSignal sig;
Vector p;
double w;
public:
prodpit() = default;
prodpit(const ProductQuadrature &q, int j0, int l);
prodpit(const prodpit &ppit) = default;
~prodpit() = default;
bool operator==(const prodpit &ppit) const;
bool
operator!=(const prodpit &ppit) const
{
return !operator==(ppit);
}
prodpit &operator=(const prodpit &spit) = delete;
prodpit &operator++();
const ParameterSignal &
signal() const
{
return sig;
}
const Vector &
point() const
{
return p;
}
double
weight() const
{
return w;
}
void print() const;
protected:
void setPointAndWeight();
};
/* The product quadrature is just QuadratureImpl with the product iterator
plugged in. The object is constructed by just giving the underlying one
dimensional quadrature, and the dimension. The only extra method is
designLevelForEvals() which for the given maximum number of evaluations (and
dimension and underlying quadrature from the object) returns a maximum level
yeilding number of evaluations less than the given number. */
class ProductQuadrature : public QuadratureImpl<prodpit>
{
friend class prodpit;
const OneDQuadrature &uquad;
public:
ProductQuadrature(int d, const OneDQuadrature &uq);
~ProductQuadrature() override = default;
int
numEvals(int l) const override
{
int res = 1;
for (int i = 0; i < dimen(); i++)
res *= uquad.numPoints(l);
return res;
}
void designLevelForEvals(int max_eval, int &lev, int &evals) const;
protected:
prodpit begin(int ti, int tn, int level) const override;
};
#endif
dynare++/integ/cc/quadrature.hh deleted 100644 → 0
View file @ 53ae5493
/*
* Copyright © 2005 Ondra Kamenik
* Copyright © 2019 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/>.
*/
// Quadrature.
/* This file defines an interface for one dimensional (non-nested) quadrature
OneDQuadrature, and a parent for all multi-dimensional quadratures. This
parent class Quadrature presents a general concept of quadrature, this is
N
∫ f(x)dx ≃ ∑ wᵢxᵢ
ⁱ⁼¹
The class Quadrature just declares this concept. The concept is implemented
by class QuadratureImpl which paralelizes the summation. All implementations
therefore wishing to use the parallel implementation should inherit from
QuadratureImpl and integration is done.
The integration concept relies on a point iterator, which goes through all
xᵢ and wᵢ for i=1,…,N. All the iterators must be able to go through only a
portion of the set i=1,…,N. This enables us to implement paralelism, for two
threads for example, one iterator goes from the beginning to the
(approximately) half, and the other goes from the half to the end.
Besides this concept of the general quadrature, this file defines also one
dimensional quadrature, which is basically a scheme of points and weights
for different levels. The class OneDQuadrature is a parent of all such
objects, the classes GaussHermite and GaussLegendre are specific
implementations for Gauss-Hermite and Gauss-Legendre quadratures resp. */
#ifndef QUADRATURE_H
#define QUADRATURE_H
#include <iostream>
#include <fstream>
#include <iomanip>
#include "vector_function.hh"
#include "int_sequence.hh"
#include "sthread.hh"
/* This pure virtual class represents a concept of one-dimensional (non-nested)
quadrature. So, one dimensional quadrature must return number of levels,
number of points in a given level, and then a point and a weight in a given
level and given order. */
class OneDQuadrature
{
public:
virtual ~OneDQuadrature() = default;
virtual int numLevels() const = 0;
virtual int numPoints(int level) const = 0;
virtual double point(int level, int i) const = 0;
virtual double weight(int lelel, int i) const = 0;
};
/* This is a general concept of multidimensional quadrature. at this general
level, we maintain only a dimension, and declare virtual functions for
integration. The function take two forms; first takes a constant
VectorFunction as an argument, creates locally VectorFunctionSet and do
calculation, second one takes as an argument VectorFunctionSet.
Part of the interface is a method returning a number of evaluations for a
specific level. Note two things: this assumes that the number of evaluations
is known apriori and thus it is not applicable for adaptive quadratures,
second for Monte Carlo type of quadrature, the level is a number of
evaluations. */
class Quadrature
{
protected:
int dim;
public:
Quadrature(int d) : dim(d)
{
}
virtual ~Quadrature() = default;
int
dimen() const
{
return dim;
}
virtual void integrate(const VectorFunction &func, int level,
int tn, Vector &out) const = 0;
virtual void integrate(VectorFunctionSet &fs, int level, Vector &out) const = 0;
virtual int numEvals(int level) const = 0;
};
/* This is just an integration worker, which works over a given QuadratureImpl.
It also needs the function, level, a specification of the subgroup of
points, and output vector.
See QuadratureImpl class declaration for details. */
template <typename _Tpit>
class QuadratureImpl;
template <typename _Tpit>
class IntegrationWorker : public sthread::detach_thread
{
const QuadratureImpl<_Tpit> &quad;
VectorFunction &func;
int level;
int ti;
int tn;
Vector &outvec;
public:
IntegrationWorker(const QuadratureImpl<_Tpit> &q, VectorFunction &f, int l,
int tii, int tnn, Vector &out)
: quad(q), func(f), level(l), ti(tii), tn(tnn), outvec(out)
{
}
/* This integrates the given portion of the integral. We obtain first and
last iterators for the portion (‘beg’ and ‘end’). Then we iterate through
the portion. and finally we add the intermediate result to the result
‘outvec’.
This method just everything up as it is coming. This might be imply large
numerical errors, perhaps in future something smarter should be implemented. */
void
operator()(std::mutex &mut) override
{
_Tpit beg = quad.begin(ti, tn, level);
_Tpit end = quad.begin(ti+1, tn, level);
Vector tmpall(outvec.length());
tmpall.zeros();
Vector tmp(outvec.length());
/* Note that since ‘beg’ came from begin(), it has empty signal and first
evaluation gets no signal */
for (_Tpit run = beg; run != end; ++run)
{
func.eval(run.point(), run.signal(), tmp);
tmpall.add(run.weight(), tmp);
}
{
std::unique_lock<std::mutex> lk{mut};
outvec.add(1.0, tmpall);
}
}
};
/* This is the class which implements the integration. The class is templated
by the iterator type. We declare a method begin() returning an iterator to
the beginnning of the ‘ti’-th portion out of total ‘tn’ portions for a given
level.
In addition, we define a method which saves all the points to a given file.
Only for debugging purposes. */
template <typename _Tpit>
class QuadratureImpl : public Quadrature
{
friend class IntegrationWorker<_Tpit>;
public:
QuadratureImpl(int d) : Quadrature(d)
{
}
/* Just fill a thread group with workes and run it. */
void
integrate(VectorFunctionSet &fs, int level, Vector &out) const override
{
// TODO: out.length()==func.outdim()
// TODO: dim == func.indim()
out.zeros();
sthread::detach_thread_group gr;
for (int ti = 0; ti < fs.getNum(); ti++)
gr.insert(std::make_unique<IntegrationWorker<_Tpit>>(*this, fs.getFunc(ti),
level, ti, fs.getNum(), out));
gr.run();
}
void
integrate(const VectorFunction &func,
int level, int tn, Vector &out) const override
{
VectorFunctionSet fs(func, tn);
integrate(fs, level, out);
}
/* Just for debugging. */
void
savePoints(const std::string &fname, int level) const
{
std::ofstream fd{fname, std::ios::out | std::ios::trunc};
if (fd.fail())
{
// TODO: raise
std::cerr << "Cannot open file " << fname << " for writing." << std::endl;
exit(EXIT_FAILURE);
}
_Tpit beg = begin(0, 1, level);
_Tpit end = begin(1, 1, level);
fd << std::setprecision(12);
for (_Tpit run = beg; run != end; ++run)
{
fd << std::setw(16) << run.weight();
for (int i = 0; i < dimen(); i++)
fd << '\t' << std::setw(16) << run.point()[i];
fd << std::endl;
}
fd.close();
}
_Tpit
start(int level) const
{
return begin(0, 1, level);
}
_Tpit
end(int level) const
{
return begin(1, 1, level);
}
protected:
virtual _Tpit begin(int ti, int tn, int level) const = 0;
};
/* This is only an interface to a precalculated data in file
precalc_quadrature.hh which is basically C coded static data. It implements
OneDQuadrature. The data file is supposed to define the following data:
number of levels, array of number of points at each level, an array of
weights and array of points. The both latter array store data level by
level. An offset for a specific level is stored in ‘offsets’ integer
sequence.
The implementing subclasses just fill the necessary data from the file, the
rest is calculated here. */
class OneDPrecalcQuadrature : public OneDQuadrature
{
int num_levels;
const int *num_points;
const double *weights;
const double *points;
IntSequence offsets;
public:
OneDPrecalcQuadrature(int nlevels, const int *npoints,
const double *wts, const double *pts)
: num_levels(nlevels), num_points(npoints),
weights(wts), points(pts), offsets(num_levels)
{
calcOffsets();
}
~OneDPrecalcQuadrature() override = default;
int
numLevels() const override
{
return num_levels;
}
int
numPoints(int level) const override
{
return num_points[level-1];
}
double
point(int level, int i) const override
{
return points[offsets[level-1]+i];
}
double
weight(int level, int i) const override
{
return weights[offsets[level-1]+i];
}
protected:
void calcOffsets();
};
/* Just precalculated Gauss-Hermite quadrature. This quadrature integrates
exactly integrals
+∞
∫ xᵏe^{−x²}dx
−∞
for level k.
Note that if pluging this one-dimensional quadrature to product or Smolyak
rule in order to integrate a function f through normally distributed inputs,
one has to wrap f to GaussConverterFunction and apply the product or Smolyak
rule to the new function.
Check precalc_quadrature.hh for available levels. */
class GaussHermite : public OneDPrecalcQuadrature
{
public:
GaussHermite();
};
/* Just precalculated Gauss-Legendre quadrature. This quadrature integrates
exactly integrals
₁
∫ xᵏdx
⁰
for level k.
Check precalc_quadrature.hh for available levels. */
class GaussLegendre : public OneDPrecalcQuadrature
{
public:
GaussLegendre();
};
#endif
dynare++/integ/cc/quasi_mcarlo.cc deleted 100644 → 0
View file @ 53ae5493
/*
* Copyright © 2005 Ondra Kamenik
* Copyright © 2019 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/>.
*/
#include "quasi_mcarlo.hh"
#include <cmath>
#include <iostream>
#include <iomanip>
#include <array>
/* Here in the constructor, we have to calculate a maximum length of ‘coeff’
array for a given ‘base’ and given maximum ‘maxn’. After allocation, we
calculate the coefficients. */
RadicalInverse::RadicalInverse(int n, int b, int mxn)
: num(n), base(b), maxn(mxn),
coeff(static_cast<int>(floor(log(static_cast<double>(maxn))/log(static_cast<double>(b)))+2), 0)
{
int nr = num;
j = -1;
do
{
j++;
coeff[j] = nr % base;
nr = nr / base;
}
while (nr > 0);
}
/* This evaluates the radical inverse. If there was no permutation, we have to
calculate:
c₀ c₁ cⱼ
── + ── + … + ────
b b² bʲ⁺¹
which is evaluated as:
⎛ ⎛⎛cⱼ 1 cⱼ₋₁⎞ 1 cⱼ₋₂⎞ ⎞ 1 c₀
⎢…⎢⎢──·─ + ────⎥·─ + ────⎥…⎥·─ + ──
⎝ ⎝⎝ b b b ⎠ b b ⎠ ⎠ b b
Now with permutation π, we have:
⎛ ⎛⎛π(cⱼ) 1 π(cⱼ₋₁)⎞ 1 π(cⱼ₋₂)⎞ ⎞ 1 π(c₀)
⎢…⎢⎢─────·─ + ───────⎥·─ + ───────⎥…⎥·─ + ─────
⎝ ⎝⎝ b b b ⎠ b b ⎠ ⎠ b b
*/
double
RadicalInverse::eval(const PermutationScheme &p) const
{
double res = 0;
for (int i = j; i >= 0; i--)
{
int cper = p.permute(i, base, coeff[i]);
res = (cper + res)/base;
}
return res;
}
/* We just add 1 to the lowest coefficient and check for overflow with respect
to the base. */
void
RadicalInverse::increase()
{
// TODO: raise if num+1 > maxn
num++;
int i = 0;
coeff[i]++;
while (coeff[i] == base)
{
coeff[i] = 0;
coeff[++i]++;
}
if (i > j)
j = i;
}
/* Debug print. */
void
RadicalInverse::print() const
{
std::cout << "n=" << num << " b=" << base << " c=";
coeff.print();
}
/* Here we have the first 170 primes. This means that we are not able to
integrate dimensions greater than 170. */
std::array<int, 170> HaltonSequence::primes = {
2, 3, 5, 7, 11, 13, 17, 19, 23, 29,
31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
73, 79, 83, 89, 97, 101, 103, 107, 109, 113,
127, 131, 137, 139, 149, 151, 157, 163, 167, 173,
179, 181, 191, 193, 197, 199, 211, 223, 227, 229,
233, 239, 241, 251, 257, 263, 269, 271, 277, 281,
283, 293, 307, 311, 313, 317, 331, 337, 347, 349,
353, 359, 367, 373, 379, 383, 389, 397, 401, 409,
419, 421, 431, 433, 439, 443, 449, 457, 461, 463,
467, 479, 487, 491, 499, 503, 509, 521, 523, 541,
547, 557, 563, 569, 571, 577, 587, 593, 599, 601,
607, 613, 617, 619, 631, 641, 643, 647, 653, 659,
661, 673, 677, 683, 691, 701, 709, 719, 727, 733,
739, 743, 751, 757, 761, 769, 773, 787, 797, 809,
811, 821, 823, 827, 829, 839, 853, 857, 859, 863,
877, 881, 883, 887, 907, 911, 919, 929, 937, 941,
947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013
};
/* This takes first ‘dim’ primes and constructs ‘dim’ radical inverses and
calls eval(). */
HaltonSequence::HaltonSequence(int n, int mxn, int dim, const PermutationScheme &p)
: num(n), maxn(mxn), per(p), pt(dim)
{
// TODO: raise if dim > num_primes
// TODO: raise if n > mxn
for (int i = 0; i < dim; i++)
ri.emplace_back(num, primes[i], maxn);
eval();
}
/* This calls RadicalInverse::increase() for all radical inverses and calls
eval(). */
void
HaltonSequence::increase()
{
for (auto & i : ri)
i.increase();
num++;
if (num <= maxn)
eval();
}
/* This sets point ‘pt’ to radical inverse evaluations in each dimension. */
void
HaltonSequence::eval()
{
for (unsigned int i = 0; i < ri.size(); i++)
pt[i] = ri[i].eval(per);
}
/* Debug print. */
void
HaltonSequence::print() const
{
auto ff = std::cout.flags();
for (const auto & i : ri)
i.print();
std::cout << "point=[ "
<< std::fixed << std::setprecision(6);
for (unsigned int i = 0; i < ri.size(); i++)
std::cout << std::setw(7) << pt[i] << ' ';
std::cout << ']' << std::endl;
std::cout.flags(ff);
}
qmcpit::qmcpit(const QMCSpecification &s, int n)
: spec(s), halton{n, s.level(), s.dimen(), s.getPerScheme()},
sig{s.dimen()}
{
}
bool
qmcpit::operator==(const qmcpit &qpit) const
{
return &spec == &qpit.spec && halton.getNum() == qpit.halton.getNum();
}
qmcpit &
qmcpit::operator++()
{
halton.increase();
return *this;
}
double
qmcpit::weight() const
{
return 1.0/spec.level();
}
int
WarnockPerScheme::permute(int i, int base, int c) const
{
return (c+i) % base;
}
int
ReversePerScheme::permute(int i, int base, int c) const
{
return (base-c) % base;
}
dynare++/integ/cc/quasi_mcarlo.hh deleted 100644 → 0
View file @ 53ae5493
/*
* Copyright © 2005 Ondra Kamenik
* Copyright © 2019 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/>.
*/
// Quasi Monte Carlo quadrature.
/* This defines quasi Monte Carlo quadratures for cube and for a function
multiplied by normal density. The quadrature for a cube is named
QMCarloCubeQuadrature and integrates:
∫ f(x)dx
[0,1]ⁿ
The quadrature for a function of normally distributed parameters is named
QMCarloNormalQuadrature and integrates:
1
──────── ∫ f(x)e^{−½xᵀx}dx
√{(2π)ⁿ} [−∞,+∞]ⁿ
For a cube we define qmcpit as iterator of QMCarloCubeQuadrature, and for
the normal density multiplied function we define qmcnpit as iterator of
QMCarloNormalQuadrature.
The quasi Monte Carlo method generates low discrepancy points with equal
weights. The one dimensional low discrepancy sequences are generated by
RadicalInverse class, the sequences are combined for higher dimensions by
HaltonSequence class. The Halton sequence can use a permutation scheme;
PermutattionScheme is an abstract class for all permutaton schemes. We have
three implementations: WarnockPerScheme, ReversePerScheme, and
IdentityPerScheme. */
#ifndef QUASI_MCARLO_H
#define QUASI_MCARLO_H
#include "int_sequence.hh"
#include "quadrature.hh"
#include "Vector.hh"
#include <vector>
/* This abstract class declares permute() method which permutes coefficient ‘c’
having index of ‘i’ fro the base ‘base’ and returns the permuted coefficient
which must be in 0,…,base−1. */
class PermutationScheme
{
public:
PermutationScheme() = default;
virtual ~PermutationScheme() = default;
virtual int permute(int i, int base, int c) const = 0;
};
/* This class represents an integer number ‘num’ as c₀+c₁b+c₂b²+…+cⱼbʲ, where b
is ‘base’ and c₀,…,cⱼ is stored in ‘coeff’. The size of IntSequence coeff
does not grow with growing ‘num’, but is fixed from the very beginning and
is set according to supplied maximum ‘maxn’.
The basic method is eval() which evaluates the RadicalInverse with a given
permutation scheme and returns the point, and increase() which increases
‘num’ and recalculates the coefficients. */
class RadicalInverse
{
int num;
int base;
int maxn;
int j;
IntSequence coeff;
public:
RadicalInverse(int n, int b, int mxn);
RadicalInverse(const RadicalInverse &ri) = default;
RadicalInverse &operator=(const RadicalInverse &radi) = default;
double eval(const PermutationScheme &p) const;
void increase();
void print() const;
};
/* This is a vector of RadicalInverses, each RadicalInverse has a different
prime as its base. The static members ‘primes’ and ‘num_primes’ define a
precalculated array of primes. The increase() method of the class increases
indices in all RadicalInverse’s and sets point ‘pt’ to contain the points in
each dimension. */
class HaltonSequence
{
private:
static std::array<int, 170> primes;
protected:
int num;
int maxn;
std::vector<RadicalInverse> ri;
const PermutationScheme &per;
Vector pt;
public:
HaltonSequence(int n, int mxn, int dim, const PermutationScheme &p);
HaltonSequence(const HaltonSequence &hs) = default;
HaltonSequence &operator=(const HaltonSequence &hs) = delete;
void increase();
const Vector &
point() const
{
return pt;
}
const int
getNum() const
{
return num;
}
void print() const;
protected:
void eval();
};
/* This is a specification of quasi Monte Carlo quadrature. It consists of
dimension ‘dim’, number of points (or level) ‘lev’, and the permutation
scheme. This class is common to all quasi Monte Carlo classes. */
class QMCSpecification
{
protected:
int dim;
int lev;
const PermutationScheme &per_scheme;
public:
QMCSpecification(int d, int l, const PermutationScheme &p)
: dim(d), lev(l), per_scheme(p)
{
}
virtual ~QMCSpecification() = default;
int
dimen() const
{
return dim;
}
int
level() const
{
return lev;
}
const PermutationScheme &
getPerScheme() const
{
return per_scheme;
}
};
/* This is an iterator for quasi Monte Carlo over a cube QMCarloCubeQuadrature.
The iterator maintains HaltonSequence of the same dimension as given by the
specification. An iterator can be constructed from a given number ‘n’, or by
a copy constructor. */
class qmcpit
{
protected:
const QMCSpecification &spec;
HaltonSequence halton;
ParameterSignal sig;
public:
qmcpit(const QMCSpecification &s, int n);
qmcpit(const qmcpit &qpit) = default;
virtual ~qmcpit() = default;
bool operator==(const qmcpit &qpit) const;
bool
operator!=(const qmcpit &qpit) const
{
return !operator==(qpit);
}
qmcpit &operator=(const qmcpit &qpit) = delete;
qmcpit &operator++();
const ParameterSignal &
signal() const
{
return sig;
}
const Vector &
point() const
{
return halton.point();
}
double weight() const;
void
print() const
{
halton.print();
}
};
/* This is an easy declaration of quasi Monte Carlo quadrature for a cube.
Everything important has been done in its iterator qmcpit, so we only
inherit from general Quadrature and reimplement begin() and numEvals(). */
class QMCarloCubeQuadrature : public QuadratureImpl<qmcpit>, public QMCSpecification
{
public:
QMCarloCubeQuadrature(int d, int l, const PermutationScheme &p)
: QuadratureImpl<qmcpit>(d), QMCSpecification(d, l, p)
{
}
~QMCarloCubeQuadrature() override = default;
int
numEvals(int l) const override
{
return l;
}
protected:
qmcpit
begin(int ti, int tn, int lev) const override
{
return qmcpit(*this, ti*level()/tn + 1);
}
};
/* Declares Warnock permutation scheme. */
class WarnockPerScheme : public PermutationScheme
{
public:
int permute(int i, int base, int c) const override;
};
/* Declares reverse permutation scheme. */
class ReversePerScheme : public PermutationScheme
{
public:
int permute(int i, int base, int c) const override;
};
/* Declares no permutation (identity) scheme. */
class IdentityPerScheme : public PermutationScheme
{
public:
int
permute(int i, int base, int c) const override
{
return c;
}
};
#endif
dynare++/integ/cc/smolyak.cc deleted 100644 → 0
View file @ 53ae5493
/*
* Copyright © 2005 Ondra Kamenik
* Copyright © 2019 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/>.
*/
#include "smolyak.hh"
#include "symmetry.hh"
#include <iostream>
#include <iomanip>
/* This constructs a beginning of ‘isum’ summand in ‘smolq’. We must be careful
here, since ‘isum’ can be past-the-end, so no reference to vectors in
‘smolq’ by ‘isum’ must be done in this case. */
smolpit::smolpit(const SmolyakQuadrature &q, unsigned int isum)
: smolq(q), isummand(isum),
jseq(q.dimen(), 0),
sig{q.dimen()},
p{q.dimen()}
{
if (isummand < q.numSummands())
setPointAndWeight();
}
bool
smolpit::operator==(const smolpit &spit) const
{
return &smolq == &spit.smolq && isummand == spit.isummand && jseq == spit.jseq;
}
/* We first try to increase index within the current summand. If we are at
maximum, we go to a subsequent summand. Note that in this case all indices
in ‘jseq’ will be zero, so no change is needed. */
smolpit &
smolpit::operator++()
{
const IntSequence &levpts = smolq.levpoints[isummand];
int i = smolq.dimen()-1;
jseq[i]++;
while (i >= 0 && jseq[i] == levpts[i])
{
jseq[i] = 0;
i--;
if (i >= 0)
jseq[i]++;
}
sig.signalAfter(std::max(i, 0));
if (i < 0)
isummand++;
if (isummand < smolq.numSummands())
setPointAndWeight();
return *this;
}
/* Here we set the point coordinates according to ‘jseq’ and
‘isummand’. Also the weight is set here. */
void
smolpit::setPointAndWeight()
{
// todo: raise if isummand ≥ smolq.numSummands()
int l = smolq.level;
int d = smolq.dimen();
int sumk = (smolq.levels[isummand]).sum();
int m1exp = l + d - sumk - 1;
w = (2*(m1exp/2) == m1exp) ? 1.0 : -1.0;
w *= PascalTriangle::noverk(d-1, sumk-l);
for (int i = 0; i < d; i++)
{
int ki = (smolq.levels[isummand])[i];
p[i] = (smolq.uquad).point(ki, jseq[i]);
w *= (smolq.uquad).weight(ki, jseq[i]);
}
}
/* Debug print. */
void
smolpit::print() const
{
auto ff = std::cout.flags();
std::cout << "isum=" << std::left << std::setw(3) << isummand << std::right << ": [";
for (int i = 0; i < smolq.dimen(); i++)
std::cout << std::setw(2) << (smolq.levels[isummand])[i] << ' ';
std::cout << "] j=[";
for (int i = 0; i < smolq.dimen(); i++)
std::cout << std::setw(2) << jseq[i] << ' ';
std::cout << std::showpos << std::fixed << std::setprecision(3)
<< "] " << std::setw(4) << w << "*(";
for (int i = 0; i < smolq.dimen()-1; i++)
std::cout << std::setw(4) << p[i] << ' ';
std::cout << std::setw(4) << p[smolq.dimen()-1] << ')' << std::endl;
std::cout.flags(ff);
}
/* Here is the constructor of SmolyakQuadrature. We have to setup ‘levels’,
‘levpoints’ and ‘cumevals’. We have to go through all d-dimensional
sequences k, such that l≤|k|≤l+d−1 and all kᵢ are positive integers. This is
equivalent to going through all k such that l−d≤|k|≤l−1 and all kᵢ are
non-negative integers. This is equivalent to going through d+1 dimensional
sequences (k,x) such that |(k,x)|=l−1 and x=0,…,d−1. The resulting sequence
of positive integers is obtained by adding 1 to all kᵢ. */
SmolyakQuadrature::SmolyakQuadrature(int d, int l, const OneDQuadrature &uq)
: QuadratureImpl<smolpit>(d), level(l), uquad(uq)
{
// TODO: check l>1, l≥d
// TODO: check l≥uquad.miLevel(), l≤uquad.maxLevel()
int cum = 0;
for (const auto &si : SymmetrySet(l-1, d+1))
{
if (si[d] <= d-1)
{
IntSequence lev(si, 0, d);
lev.add(1);
levels.push_back(lev);
IntSequence levpts(d);
for (int i = 0; i < d; i++)
levpts[i] = uquad.numPoints(lev[i]);
levpoints.push_back(levpts);
cum += levpts.mult();
cumevals.push_back(cum);
}
}
}
/* Here we return a number of evalutions of the quadrature for the given level.
If the given level is the current one, we simply return the maximum
cumulative number of evaluations. Otherwise we call costly
calcNumEvaluations() method. */
int
SmolyakQuadrature::numEvals(int l) const
{
if (l != level)
return calcNumEvaluations(l);
else
return cumevals[numSummands()-1];
}
/* This divides all the evaluations to ‘tn’ approximately equal groups, and
returns the beginning of the specified group ‘ti’. The granularity of
divisions are summands as listed by ‘levels’. */
smolpit
SmolyakQuadrature::begin(int ti, int tn, int l) const
{
// TODO: raise is level≠l
if (ti == tn)
return smolpit(*this, numSummands());
int totevals = cumevals[numSummands()-1];
int evals = (totevals*ti)/tn;
unsigned int isum = 0;
while (isum+1 < numSummands() && cumevals[isum+1] < evals)
isum++;
return smolpit(*this, isum);
}
/* This is the same in a structure as SmolyakQuadrature constructor. We have to
go through all summands and calculate a number of evaluations in each
summand. */
int
SmolyakQuadrature::calcNumEvaluations(int lev) const
{
int cum = 0;
for (const auto &si : SymmetrySet(lev-1, dim+1))
{
if (si[dim] <= dim-1)
{
IntSequence lev(si, 0, dim);
lev.add(1);
IntSequence levpts(dim);
for (int i = 0; i < dim; i++)
levpts[i] = uquad.numPoints(lev[i]);
cum += levpts.mult();
}
}
return cum;
}
/* This returns a maximum level such that the number of evaluations is less
than the given number. */
void
SmolyakQuadrature::designLevelForEvals(int max_evals, int &lev, int &evals) const
{
int last_evals;
evals = 1;
lev = 1;
do
{
lev++;
last_evals = evals;
evals = calcNumEvaluations(lev);
}
while (lev < uquad.numLevels() && evals <= max_evals);
lev--;
evals = last_evals;
}
dynare++/integ/cc/smolyak.hh deleted 100644 → 0
View file @ 53ae5493
/*
* Copyright © 2005 Ondra Kamenik
* Copyright © 2019 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/>.
*/
// Smolyak quadrature.
/* This file defines Smolyak (sparse grid) multidimensional quadrature for
non-nested underlying one dimensional quadrature. Let Q¹ₗ denote the one
dimensional quadrature of l level. Let nₗ denote a number of points in the l
level. Than the Smolyak quadrature can be defined as
⎛ d−1 ⎞
Qᵈf = ∑ (−1)^{l+d−|k|−1} ⎢ ⎥ (Q¹_k₁⊗…⊗Q¹_{k_d})f
l≤|k|≤l+d−1 ⎝|k|−l⎠
where d is the dimension, k is d-dimensional sequence of integers, and |k|
denotes the sum of the sequence.
Here we define smolpit as Smolyak iterator and SmolyakQuadrature. */
#ifndef SMOLYAK_H
#define SMOLYAK_H
#include "int_sequence.hh"
#include "tl_static.hh"
#include "vector_function.hh"
#include "quadrature.hh"
#include "pascal_triangle.hh"
/* Here we define the Smolyak point iterator. The Smolyak formula can be broken
to a sum of product quadratures with various combinations of levels. The
iterator follows this pattern. It maintains an index to a summand and then a
point coordinates within the summand (product quadrature). The array of
summands to which the ‘isummand’ points is maintained by the
SmolyakQuadrature class to which the object knows the pointer ‘smolq’.
We provide a constructor which points to the beginning of the given summand.
This constructor is used in SmolyakQuadrature::begin() method which
approximately divideds all the iterators to subsets of equal size. */
class SmolyakQuadrature;
class smolpit
{
protected:
const SmolyakQuadrature &smolq;
unsigned int isummand{0};
IntSequence jseq;
ParameterSignal sig;
Vector p;
double w;
public:
smolpit(const SmolyakQuadrature &q, unsigned int isum);
smolpit(const smolpit &spit) = default;
~smolpit() = default;
bool operator==(const smolpit &spit) const;
bool
operator!=(const smolpit &spit) const
{
return !operator==(spit);
}
smolpit &operator=(const smolpit &spit) = delete;
smolpit &operator++();
const ParameterSignal &
signal() const
{
return sig;
}
const Vector &
point() const
{
return p;
}
double
weight() const
{
return w;
}
void print() const;
protected:
void setPointAndWeight();
};
/* Here we define the class SmolyakQuadrature. It maintains an array of
summands of the Smolyak quadrature formula:
⎛ d−1 ⎞
∑ (−1)^{l+d−|k|−1} ⎢ ⎥ (Q¹_k₁⊗…⊗Q¹_{k_d})f
l≤|k|≤l+d−1 ⎝|k|−l⎠
Each summand is fully specified by sequence k. The summands are here
represented (besides k) also by sequence of number of points in each level
selected by k, and also by a cummulative number of evaluations. The latter
two are added only for conveniency.
The summands in the code are given by ‘levels’, which is a vector of
k sequences, further by ‘levpoints’ which is a vector of sequences
of nuber of points in each level, and by ‘cumevals’ which is the
cumulative number of points, this is:
d
∑ ∏ n_kᵢ
ᵏ ⁱ⁼¹
where the sum is done through all k before the current.
The ‘levels’ and ‘levpoints’ vectors are used by smolpit. */
class SmolyakQuadrature : public QuadratureImpl<smolpit>
{
friend class smolpit;
int level;
const OneDQuadrature &uquad;
std::vector<IntSequence> levels;
std::vector<IntSequence> levpoints;
std::vector<int> cumevals;
public:
SmolyakQuadrature(int d, int l, const OneDQuadrature &uq);
~SmolyakQuadrature() override = default;
int numEvals(int level) const override;
void designLevelForEvals(int max_eval, int &lev, int &evals) const;
protected:
smolpit begin(int ti, int tn, int level) const override;
unsigned int
numSummands() const
{
return levels.size();
}
private:
int calcNumEvaluations(int level) const;
};
#endif
dynare++/integ/cc/vector_function.cc deleted 100644 → 0
View file @ 53ae5493
/*
* Copyright © 2005 Ondra Kamenik
* Copyright © 2019 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/>.
*/
#include "vector_function.hh"
#include <dynlapack.h>
#include <cmath>
#include <algorithm>
/* Just an easy constructor of sequence of booleans defaulting to change
everywhere. */
ParameterSignal::ParameterSignal(int n)
: data(n, true)
{
}
/* This sets ‘false’ (no change) before a given parameter, and ‘true’ (change)
after the given parameter (including). */
void
ParameterSignal::signalAfter(int l)
{
for (size_t i = 0; i < std::min(static_cast<size_t>(l), data.size()); i++)
data[i] = false;
for (size_t i = l; i < data.size(); i++)
data[i] = true;
}
/* This constructs a function set hardcopying also the first. */
VectorFunctionSet::VectorFunctionSet(const VectorFunction &f, int n)
: funcs(n)
{
for (int i = 0; i < n; i++)
{
func_copies.push_back(f.clone());
funcs[i] = func_copies.back().get();
}
}
/* This constructs a function set with shallow copy in the first and hard
copies in others. */
VectorFunctionSet::VectorFunctionSet(VectorFunction &f, int n)
: funcs(n)
{
if (n > 0)
funcs[0] = &f;
for (int i = 1; i < n; i++)
{
func_copies.push_back(f.clone());
funcs[i] = func_copies.back().get();
}
}
/* Here we construct the object from the given function f and given
variance-covariance matrix Σ=vcov. The matrix A is calculated as lower
triangular and yields Σ=AAᵀ. */
GaussConverterFunction::GaussConverterFunction(VectorFunction &f, const GeneralMatrix &vcov)
: VectorFunction(f), func(&f), A(vcov.nrows(), vcov.nrows()),
multiplier(calcMultiplier())
{
// TODO: raise if A.nrows() ≠ indim()
calcCholeskyFactor(vcov);
}
GaussConverterFunction::GaussConverterFunction(std::unique_ptr<VectorFunction> f, const GeneralMatrix &vcov)
: VectorFunction(*f), func_storage{move(f)}, func{func_storage.get()}, A(vcov.nrows(), vcov.nrows()),
multiplier(calcMultiplier())
{
// TODO: raise if A.nrows() ≠ indim()
calcCholeskyFactor(vcov);
}
GaussConverterFunction::GaussConverterFunction(const GaussConverterFunction &f)
: VectorFunction(f), func_storage{f.func->clone()}, func{func_storage.get()}, A(f.A),
multiplier(f.multiplier)
{
}
/* Here we evaluate the function
g(y) = 1/√(πⁿ) f(√2·Ay).
Since the matrix A is lower triangular, the change signal for the function f
will look like (0,…,0,1,…,1) where the first 1 is in the same position as
the first change in the given signal ‘sig’ of the input y=point. */
void
GaussConverterFunction::eval(const Vector &point, const ParameterSignal &sig, Vector &out)
{
ParameterSignal s(sig);
int i = 0;
while (i < indim() && !sig[i])
i++;
s.signalAfter(i);
Vector x(indim());
x.zeros();
A.multaVec(x, point);
x.mult(sqrt(2.0));
func->eval(x, s, out);
out.mult(multiplier);
}
/* This returns 1/√(πⁿ). */
double
GaussConverterFunction::calcMultiplier() const
{
return sqrt(pow(M_PI, -1*indim()));
}
void
GaussConverterFunction::calcCholeskyFactor(const GeneralMatrix &vcov)
{
A = vcov;
lapack_int rows = A.nrows(), lda = A.getLD();
for (int i = 0; i < rows; i++)
for (int j = i+1; j < rows; j++)
A.get(i, j) = 0.0;
lapack_int info;
dpotrf("L", &rows, A.base(), &lda, &info);
// TODO: raise if info≠1
}
dynare++/integ/cc/vector_function.hh deleted 100644 → 0
View file @ 53ae5493
/*
* Copyright © 2005 Ondra Kamenik
* Copyright © 2019 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/>.
*/
// Vector function.
/* This file defines interface for functions taking a vector as an input and
returning a vector (with a different size) as an output. We are also
introducing a parameter signalling; it is a boolean vector which tracks
parameters which were changed from the previous call. The VectorFunction
implementation can exploit this information and evaluate the function more
efficiently. The information can be completely ignored.
From the signalling reason, and from other reasons, the function evaluation
is not const. */
#ifndef VECTOR_FUNCTION_H
#define VECTOR_FUNCTION_H
#include "Vector.hh"
#include "GeneralMatrix.hh"
#include <vector>
#include <memory>
/* This is a simple class representing a vector of booleans. The items night be
retrieved or changed, or can be set ‘true’ after some point. This is useful
when we multiply the vector with lower triangular matrix.
‘true’ means that a parameter was changed. */
class ParameterSignal
{
protected:
std::vector<bool> data;
public:
ParameterSignal(int n);
ParameterSignal(const ParameterSignal &sig) = default;
~ParameterSignal() = default;
void signalAfter(int l);
bool
operator[](int i) const
{
return data[i];
}
std::vector<bool>::reference
operator[](int i)
{
return data[i];
}
};
/* This is the abstract class for vector function. At this level of abstraction
we only need to know size of input vector and a size of output vector.
The important thing here is a clone method, we will need to make hard copies
of vector functions since the evaluations are not const. The hardcopies
apply for parallelization. */
class VectorFunction
{
protected:
int in_dim;
int out_dim;
public:
VectorFunction(int idim, int odim)
: in_dim(idim), out_dim(odim)
{
}
VectorFunction(const VectorFunction &func) = default;
virtual ~VectorFunction() = default;
virtual std::unique_ptr<VectorFunction> clone() const = 0;
virtual void eval(const Vector &point, const ParameterSignal &sig, Vector &out) = 0;
int
indim() const
{
return in_dim;
}
int
outdim() const
{
return out_dim;
}
};
/* This makes ‘n’ copies of VectorFunction. The first constructor make exactly
‘n’ new copies, the second constructor copies only the pointer to the first
and others are hard (real) copies.
The class is useful for making a given number of copies at once, and this
set can be reused many times if we need mupliple copis of the function (for
example for paralelizing the code). */
class VectorFunctionSet
{
private:
// Stores the hard copies made by the class
std::vector<std::unique_ptr<VectorFunction>> func_copies;
protected:
std::vector<VectorFunction *> funcs;
public:
VectorFunctionSet(const VectorFunction &f, int n);
VectorFunctionSet(VectorFunction &f, int n);
~VectorFunctionSet() = default;
VectorFunction &
getFunc(int i)
{
return *(funcs[i]);
}
int
getNum() const
{
return funcs.size();
}
};
/* This class wraps another VectorFunction to allow integration of a function
through normally distributed inputs. Namely, if one wants to integrate
1
─────────── ∫ f(x)e^{−½xᵀ|Σ|⁻¹x}dx
√{(2π)ⁿ|Σ|}
then if we write Σ=AAᵀ and x=√2·Ay, we get integral
1 1
─────────── ∫ f(√2·Ay)e^{−½yᵀy} √(2ⁿ)|A|dy = ───── ∫ f(√2·Ay)e^{−½yᵀy}dy
√{(2π)ⁿ|Σ|} √(πⁿ)
which means that a given function f we have to wrap to yield a function
g(y) = 1/√(πⁿ) f(√2·Ay).
This is exactly what this class is doing. This transformation is useful
since the Gauss-Hermite points and weights are defined for weighting
function e^{−y²}, so this transformation allows using Gauss-Hermite
quadratures seemlessly in a context of integration through normally
distributed inputs.
The class maintains a pointer to the function f. When the object is
constructed by the first constructor, the f is assumed to be owned by the
caller. If the object of this class is copied, then f is copied and hence
stored in a std::unique_ptr. The second constructor takes a smart pointer to
the function and in that case the class takes ownership of f. */
class GaussConverterFunction : public VectorFunction
{
private:
std::unique_ptr<VectorFunction> func_storage;
protected:
VectorFunction *func;
GeneralMatrix A;
double multiplier;
public:
GaussConverterFunction(VectorFunction &f, const GeneralMatrix &vcov);
GaussConverterFunction(std::unique_ptr<VectorFunction> f, const GeneralMatrix &vcov);
GaussConverterFunction(const GaussConverterFunction &f);
~GaussConverterFunction() override = default;
std::unique_ptr<VectorFunction>
clone() const override
{
return std::make_unique<GaussConverterFunction>(*this);
}
void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override;
private:
double calcMultiplier() const;
void calcCholeskyFactor(const GeneralMatrix &vcov);
};
#endif
dynare++/integ/src/Makefile.am deleted 100644 → 0
View file @ 53ae5493
noinst_PROGRAMS = quadrature-points
quadrature_points_SOURCES = quadrature-points.cc
quadrature_points_CPPFLAGS = -I../.. -I../../sylv/cc -I../../integ/cc -I../../tl/cc -I../../utils/cc
quadrature_points_CXXFLAGS = $(AM_CXXFLAGS) $(THREAD_CXXFLAGS)
quadrature_points_LDADD = ../cc/libinteg.a ../../tl/cc/libtl.a ../../parser/cc/libparser.a ../../sylv/cc/libsylv.a ../../utils/cc/libutils.a $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS)