QuantEcon

Represents a scalar ARMA(p, q) process

If $\phi$ and $\theta$ are scalars, then the model is understood to be

\[ X_t = \phi X_{t-1} + \epsilon_t + \theta \epsilon_{t-1}\]

where $\epsilon_t$ is a white noise process with standard deviation sigma.

If $\phi$ and $\theta$ are arrays or sequences, then the interpretation is the ARMA(p, q) model

\[ X_t = \phi_1 X_{t-1} + ... + \phi_p X_{t-p} + \epsilon_t + \theta_1 \epsilon_{t-1} + \ldots + \theta_q \epsilon_{t-q}\]

where

• $\phi = (\phi_1, \phi_2, \ldots , \phi_p)$
• $\theta = (\theta_1, \theta_2, \ldots , \theta_q)$
• $\sigma$ is a scalar, the standard deviation of the white noise

Fields

• phi::Vector : AR parameters $\phi_1, \ldots, \phi_p$
• theta::Vector : MA parameters $\theta_1, \ldots, \theta_q$
• p::Integer : Number of AR coefficients
• q::Integer : Number of MA coefficients
• sigma::Real : Standard deviation of white noise
• ma_poly::Vector : MA polynomial –- filtering representatoin
• ar_poly::Vector : AR polynomial –- filtering representation

Examples

using QuantEcon
phi = 0.5
theta = [0.0, -0.8]
sigma = 1.0
lp = ARMA(phi, theta, sigma)
require(joinpath(dirname(@__FILE__),"..", "examples", "arma_plots.jl"))
quad_plot(lp)

Type used to evaluate constant Frisch elasticity (CFE) utility. CFE utility takes the form

v(l) = ξ l^(1 + 1/ϕ) / (1 + 1/ϕ)

Additionally, this code assumes that if l < 1e-10 then

v(l) = ξ (1e-10^(1 + 1/ϕ) / (1 + 1/ϕ) - 1e-10^(1/ϕ) * (1e-10 - l))

Type used to evaluate CRRA utility. CRRA utility takes the form

u(c) = ξ c^(1 - γ) / (1 - γ)

Additionally, this code assumes that if c < 1e-10 then

u(c) = ξ (1e-10^(1 - γ) / (1 - γ) + 1e-10^(-γ) * (c - 1e-10))

DiscreteDP type for specifying paramters for discrete dynamic programming model

Parameters

• R::Array{T,NR} : Reward Array
• Q::Array{T,NQ} : Transition Probability Array
• beta::Float64 : Discount Factor
• a_indices::Vector{Tind}: Action Indices. Empty unless using SA formulation
• a_indptr::Vector{Tind}: Action Index Pointers. Empty unless using SA formulation

Returns

• ddp::DiscreteDP : DiscreteDP object

DiscreteDP type for specifying parameters for discrete dynamic programming model Dense Matrix Formulation

Parameters

• R::Array{T,NR} : Reward Array
• Q::Array{T,NQ} : Transition Probability Array
• beta::Float64 : Discount Factor

Returns

• ddp::DiscreteDP : Constructor for DiscreteDP object

DiscreteDP type for specifying parameters for discrete dynamic programming model State-Action Pair Formulation

Parameters

• R::Array{T,NR} : Reward Array
• Q::Array{T,NQ} : Transition Probability Array
• beta::Float64 : Discount Factor
• s_indices::Vector{Tind}: State Indices. Empty unless using SA formulation
• a_indices::Vector{Tind}: Action Indices. Empty unless using SA formulation
• a_indptr::Vector{Tind}: Action Index Pointers. Empty unless using SA formulation

Returns

• ddp::DiscreteDP : Constructor for DiscreteDP object

Generates an array of draws from a discrete random variable with vector of probabilities given by q.

Fields

• q::AbstractVector: A vector of non-negative probabilities that sum to 1
• Q::AbstractVector: The cumulative sum of q

Type used to evaluate elliptical utility function. Elliptical utility takes form

v(l) = b (1 - l^μ)^(1 / μ)

A look ahead estimator associated with a given stochastic kernel p and a vector of observations X.

Fields

• p::Function: The stochastic kernel. Signature is p(x, y) and it should be vectorized in both inputs
• X::Matrix: A vector containing observations. Note that this can be passed as any kind of AbstractArray and will be coerced into an n x 1 vector.

Linear quadratic optimal control of either infinite or finite horizon

The infinite horizon problem can be written

\[\min \mathbb{E} \sum_{t=0}^{\infty} \beta^t r(x_t, u_t)\]

with

\[r(x_t, u_t) := x_t' R x_t + u_t' Q u_t + 2 u_t' N x_t\]

The finite horizon form is

\[\min \mathbb{E} \sum_{t=0}^{T-1} \beta^t r(x_t, u_t) + \beta^T x_T' R_f x_T\]

Both are minimized subject to the law of motion

\[x_{t+1} = A x_t + B u_t + C w_{t+1}\]

Here $x$ is n x 1, $u$ is k x 1, $w$ is j x 1 and the matrices are conformable for these dimensions. The sequence ${w_t}$ is assumed to be white noise, with zero mean and $\mathbb{E} w_t w_t' = I$, the j x j identity.

For this model, the time $t$ value (i.e., cost-to-go) function $V_t$ takes the form

\[x' P_T x + d_T\]

and the optimal policy is of the form $u_T = -F_T x_T$. In the infinite horizon case, $V, P, d$ and $F$ are all stationary.

Fields

• Q::ScalarOrArray : k x k payoff coefficient for control variable u. Must be symmetric and nonnegative definite
• R::ScalarOrArray : n x n payoff coefficient matrix for state variable x. Must be symmetric and nonnegative definite
• A::ScalarOrArray : n x n coefficient on state in state transition
• B::ScalarOrArray : n x k coefficient on control in state transition
• C::ScalarOrArray : n x j coefficient on random shock in state transition
• N::ScalarOrArray : k x n cross product in payoff equation
• bet::Real : Discount factor in [0, 1]
• capT::Union{Int, Void} : Terminal period in finite horizon problem
• rf::ScalarOrArray : n x n terminal payoff in finite horizon problem. Must be symmetric and nonnegative definite
• P::ScalarOrArray : n x n matrix in value function representation $V(x) = x'Px + d$
• d::Real : Constant in value function representation
• F::ScalarOrArray : Policy rule that specifies optimal control in each period

Main constructor for LQ type

Specifies default argumets for all fields not part of the payoff function or transition equation.

Arguments

• Q::ScalarOrArray : k x k payoff coefficient for control variable u. Must be symmetric and nonnegative definite
• R::ScalarOrArray : n x n payoff coefficient matrix for state variable x. Must be symmetric and nonnegative definite
• A::ScalarOrArray : n x n coefficient on state in state transition
• B::ScalarOrArray : n x k coefficient on control in state transition
• ;C::ScalarOrArray{zero(size(R}(1))) : n x j coefficient on random shock in state transition
• ;N::ScalarOrArray{zero(size(B,1)}(size(A, 2))) : k x n cross product in payoff equation
• ;bet::Real(1.0) : Discount factor in [0, 1]
• capT::Union{Int, Void}(Void) : Terminal period in finite horizon problem
• rf::ScalarOrArray{fill(NaN}(size(R)...)) : n x n terminal payoff in finite horizon problem. Must be symmetric and nonnegative definite.

A type that describes the Gaussian Linear State Space Model of the form:

\[ x_{t+1} = A x_t + C w_{t+1} \\ y_t = G x_t + H v_t\]

where ${w_t}$ and ${v_t}$ are independent and standard normal with dimensions k and l respectively. The initial conditions are $\mu_0$ and $\Sigma_0$ for $x_0 \sim N(\mu_0, \Sigma_0)$. When $\Sigma_0=0$, the draw of $x_0$ is exactly $\mu_0$.

Fields

• A::Matrix Part of the state transition equation. It should be n x n
• C::Matrix Part of the state transition equation. It should be n x m
• G::Matrix Part of the observation equation. It should be k x n
• H::Matrix Part of the observation equation. It should be k x l
• k::Int Dimension
• n::Int Dimension
• m::Int Dimension
• l::Int Dimension
• mu_0::Vector This is the mean of initial draw and is of length n
• Sigma_0::Matrix This is the variance of the initial draw and is n x n and also should be positive definite and symmetric

Linear interpolation in one dimension

Fields

• breaks::AbstractVector : A sorted array of grid points on which to interpolate

• vals::AbstractVector : The function values associated with each of the grid points

Examples

breaks = cumsum(0.1 .* rand(20))
vals = 0.1 .* sin.(breaks)
li = LinInterp(breaks, vals)

# do interpolation via `call` method on a LinInterp object
li(0.2)

# use broadcasting to evaluate at multiple points
li.([0.1, 0.2, 0.3])

Type used to evaluate log utility. Log utility takes the form

u(c) = \log(c)

Additionally, this code assumes that if c < 1e-10 then

u(c) = log(1e-10) + 1e10*(c - 1e-10)

This refers to the Modified Policy Iteration solution algorithm.

References

https://lectures.quantecon.org/jl/discrete_dp.html

Finite-state discrete-time Markov chain.

Methods are available that provide useful information such as the stationary distributions, and communication and recurrent classes, and allow simulation of state transitions.

Fields

• p::AbstractMatrix : The transition matrix. Must be square, all elements must be nonnegative, and all rows must sum to unity.
• state_values::AbstractVector : Vector containing the values associated with the states.

Returns the controlled Markov chain for a given policy sigma.

Parameters

• ddp::DiscreteDP : Object that contains the model parameters
• ddpr::DPSolveResult : Object that contains result variables

Returns

mc : MarkovChain Controlled Markov chain.

This refers to the Policy Iteration solution algorithm.

References

https://lectures.quantecon.org/jl/discrete_dp.html

Represents infinite horizon robust LQ control problems of the form

\[ \min_{u_t} \sum_t \beta^t {x_t' R x_t + u_t' Q u_t }\]

subject to

\[ x_{t+1} = A x_t + B u_t + C w_{t+1}\]

and with model misspecification parameter $\theta$.

Fields

• Q::Matrix{Float64} : The cost(payoff) matrix for the controls. See above for more. $Q$ should be k x k and symmetric and positive definite
• R::Matrix{Float64} : The cost(payoff) matrix for the state. See above for more. $R$ should be n x n and symmetric and non-negative definite
• A::Matrix{Float64} : The matrix that corresponds with the state in the state space system. $A$ should be n x n
• B::Matrix{Float64} : The matrix that corresponds with the control in the state space system. $B$ should be n x k
• C::Matrix{Float64} : The matrix that corresponds with the random process in the state space system. $C$ should be n x j
• beta::Real : The discount factor in the robust control problem
• theta::Real The robustness factor in the robust control problem
• k, n, j::Int : Dimensions of input matrices

SimplexGrid

Iterator version of simplex_grid, i.e., iterator that iterates over the integer points in the (m-1)-dimensional simplex $\{x \mid x_1 + \cdots + x_m = n, x_i \geq 0\}$, or equivalently, the m-part compositions of n, in lexicographic order.

Fields

• m::Int : Dimension of each point. Must be a positive integer.
• n::Int : Number which the coordinates of each point sum to. Must be a nonnegative integer.

Examples

julia> sg = SimplexGrid(3, 4);

julia> for x in sg
           @show x
       end
x = [0, 0, 4]
x = [0, 1, 3]
x = [0, 2, 2]
x = [0, 3, 1]
x = [0, 4, 0]
x = [1, 0, 3]
x = [1, 1, 2]
x = [1, 2, 1]
x = [1, 3, 0]
x = [2, 0, 2]
x = [2, 1, 1]
x = [2, 2, 0]
x = [3, 0, 1]
x = [3, 1, 0]
x = [4, 0, 0]

This refers to the Value Iteration solution algorithm.

References

https://lectures.quantecon.org/jl/discrete_dp.html

Computes the periodogram

\[I(w) = \frac{1}{n} | \sum_{t=0}^{n-1} x_t e^{itw} |^2\]

at the Fourier frequences $w_j := 2 \frac{\pi j}{n}, j = 0, \ldots, n - 1$, using the fast Fourier transform. Only the frequences $w_j$ in $[0, \pi]$ and corresponding values $I(w_j)$ are returned. If a window type is given then smoothing is performed.

Arguments

• x::Array: An array containing the data to smooth
• window_len::Int(7): An odd integer giving the length of the window
• window::AbstractString("hanning"): A string giving the window type. Possible values are flat, hanning, hamming, bartlett, or blackman

Returns

• w::Array{Float64}: Fourier frequencies at which the periodogram is evaluated
• I_w::Array{Float64}: The periodogram at frequences w

Return the period of the Markov chain mc.

Arguments

• mc::MarkovChain : MarkovChain instance.

Returns

• ::Int : Period of mc.

Compute agent 2's best cost-minimizing response $K$, given $F$.

Arguments

• rlq::RBLQ: Instance of RBLQ type
• F::Matrix{Float64}: A k x n array representing agent 1's policy

Returns

• K::Matrix{Float64} : Agent's best cost minimizing response corresponding to $F$
• P::Matrix{Float64} : The value function corresponding to $F$

Compute agent 1's best cost-minimizing response $K$, given $F$.

Arguments

• rlq::RBLQ: Instance of RBLQ type
• K::Matrix{Float64}: A k x n array representing the worst case matrix

Returns

• F::Matrix{Float64} : Agent's best cost minimizing response corresponding to $K$
• P::Matrix{Float64} : The value function corresponding to $K$

Method of RQ_sigma that extracts sigma from a DPSolveResult

See other docstring for details

Given a policy sigma, return the reward vector R_sigma and the transition probability matrix Q_sigma.

Parameters

• ddp::DiscreteDP : Object that contains the model parameters
• sigma::AbstractVector{Int}: policy rule vector

Returns

• R_sigma::Array{Float64}: Reward vector for sigma, of length n.

• Q_sigma::Array{Float64}: Transition probability matrix for sigma, of shape (n, n).

Compute periodogram from data x, using prewhitening, smoothing and recoloring. The data is fitted to an AR(1) model for prewhitening, and the residuals are used to compute a first-pass periodogram with smoothing. The fitted coefficients are then used for recoloring.

Arguments

• x::Array: An array containing the data to smooth
• window_len::Int(7): An odd integer giving the length of the window
• window::AbstractString("hanning"): A string giving the window type. Possible values are flat, hanning, hamming, bartlett, or blackman

Returns

• w::Array{Float64}: Fourier frequencies at which the periodogram is evaluated
• I_w::Array{Float64}: The periodogram at frequences w

Compute the autocovariance function from the ARMA parameters over the integers range(num_autocov) using the spectral density and the inverse Fourier transform.

Arguments

• arma::ARMA: Instance of ARMA type
• ;num_autocov::Integer(16) : The number of autocovariances to calculate

The $D$ operator, mapping $P$ into

\[ B(P) := R - \beta^2 A'PB(Q + \beta B'PB)^{-1}B'PA + \beta A'PA\]

and also returning

\[ F := (Q + \beta B'PB)^{-1} \beta B'PA\]

Arguments

• rlq::RBLQ: Instance of RBLQ type
• P::Matrix{Float64} : size is n x n

Returns

• F::Matrix{Float64} : The $F$ matrix as defined above
• new_p::Matrix{Float64} : The matrix $P$ after applying the $B$ operator

backward_induction(ddp, J[, v_term=zeros(num_states(ddp))])

Solve by backward induction a $J$-period finite horizon discrete dynamic program with stationary reward $r$ and transition probability functions $q$ and discount factor $\beta \in [0, 1]$.

The optimal value functions $v^{\ast}_1, \ldots, v^{\ast}_{J+1}$ and policy functions $\sigma^{\ast}_1, \ldots, \sigma^{\ast}_J$ are obtained by $v^{\ast}_{J+1} = v_{J+1}$, and

\[v^{\ast}_j(s) = \max_{a \in A(s)} r(s, a) + \beta \sum_{s' \in S} q(s'|s, a) v^{\ast}_{j+1}(s') \quad (s \in S)\]

and

\[\sigma^{\ast}_j(s) \in \operatorname*{arg\,max}_{a \in A(s)} r(s, a) + \beta \sum_{s' \in S} q(s'|s, a) v^*_{j+1}(s') \quad (s \in S)\]

for $j= J, \ldots, 1$, where the terminal value function $v_{J+1}$ is exogenously given by v_term.

Parameters

• ddp::DiscreteDP{T} : Object that contains the Model Parameters
• J::Integer: Number of decision periods
• v_term::AbstractVector{<:Real}=zeros(num_states(ddp)): Terminal value function of length equal to n (the number of states)

Returns

• vs::Matrix{S}: Array of shape (n, J+1) where vs[:,j] contains the optimal value function at period j = 1, ..., J+1.
• sigmas::Matrix{Int}: Array of shape (n, J) where sigmas[:,j] contains the optimal policy function at period j = 1, ..., J.

The Bellman operator, which computes and returns the updated value function $Tv$ for a value function $v$.

Parameters

• ddp::DiscreteDP : Object that contains the model parameters
• v::AbstractVector{T<:AbstractFloat}: The current guess of the value function
• Tv::AbstractVector{T<:AbstractFloat}: A buffer array to hold the updated value function. Initial value not used and will be overwritten
• sigma::AbstractVector: A buffer array to hold the policy function. Initial values not used and will be overwritten

Returns

• Tv::typeof(Tv) : Updated value function vector
• sigma::typeof(sigma) : Updated policiy function vector

Apply the Bellman operator using v=ddpr.v, Tv=ddpr.Tv, and sigma=ddpr.sigma

Notes

Updates ddpr.Tv and ddpr.sigma inplace

The Bellman operator, which computes and returns the updated value function $Tv$ for a given value function $v$.

This function will fill the input v with Tv and the input sigma with the corresponding policy rule.

Parameters

• ddp::DiscreteDP: The ddp model
• v::AbstractVector{T<:AbstractFloat}: The current guess of the value function. This array will be overwritten
• sigma::AbstractVector: A buffer array to hold the policy function. Initial values not used and will be overwritten

Returns

• Tv::Vector: Updated value function vector
• sigma::typeof(sigma): Policy rule

The Bellman operator, which computes and returns the updated value function $Tv$ for a given value function $v$.

Parameters

• ddp::DiscreteDP: The ddp model
• v::AbstractVector: The current guess of the value function

Returns

• Tv::Vector : Updated value function vector

Find the root of the f on the bracketing inverval [x1, x2] via bisection.

Arguments

• f::Function: The function you want to bracket
• x1::T: Lower border for search interval
• x2::T: Upper border for search interval
• ;maxiter::Int(500): Maximum number of bisection iterations
• ;xtol::Float64(1e-12): The routine converges when a root is known to lie within xtol of the value return. Should be >= 0. The routine modifies this to take into account the relative precision of doubles.
• ;rtol::Float64(2*eps()):The routine converges when a root is known to lie within rtol times the value returned of the value returned. Should be ≥ 0

Returns

• x::T: The found root

Exceptions

• Throws an ArgumentError if [x1, x2] does not form a bracketing interval
• Throws a ConvergenceError if the maximum number of iterations is exceeded

References

Matches bisect function from scipy/scipy/optimize/Zeros/bisect.c

Find the root of the f on the bracketing inverval [x1, x2] via brent's algo.

Arguments

• f::Function: The function you want to bracket
• x1::T: Lower border for search interval
• x2::T: Upper border for search interval
• ;maxiter::Int(500): Maximum number of bisection iterations
• ;xtol::Float64(1e-12): The routine converges when a root is known to lie within xtol of the value return. Should be >= 0. The routine modifies this to take into account the relative precision of doubles.
• ;rtol::Float64(2*eps()):The routine converges when a root is known to lie within rtol times the value returned of the value returned. Should be ≥ 0

Returns

• x::T: The found root

Exceptions

• Throws an ArgumentError if [x1, x2] does not form a bracketing interval
• Throws a ConvergenceError if the maximum number of iterations is exceeded

References

Matches brentq function from scipy/scipy/optimize/Zeros/bisectq.c

Find a root of the f on the bracketing inverval [x1, x2] via modified brent

This routine uses a hyperbolic extrapolation formula instead of the standard inverse quadratic formula. Otherwise it is the original Brent's algorithm, as implemented in the brent function.

Arguments

• f::Function: The function you want to bracket
• x1::T: Lower border for search interval
• x2::T: Upper border for search interval
• ;maxiter::Int(500): Maximum number of bisection iterations
• ;xtol::Float64(1e-12): The routine converges when a root is known to lie within xtol of the value return. Should be >= 0. The routine modifies this to take into account the relative precision of doubles.
• ;rtol::Float64(2*eps()):The routine converges when a root is known to lie within rtol times the value returned of the value returned. Should be ≥ 0

Returns

• x::T: The found root

Exceptions

• Throws an ArgumentError if [x1, x2] does not form a bracketing interval
• Throws a ConvergenceError if the maximum number of iterations is exceeded

References

Matches brenth function from scipy/scipy/optimize/Zeros/bisecth.c

ckron(arrays::AbstractArray...)

Repeatedly apply kronecker products to the arrays. Equilvalent to reduce(kron, arrays)

Find the communication classes of the Markov chain mc.

Arguments

• mc::MarkovChain : MarkovChain instance.

Returns

• ::Vector{Vector{Int}} : Vector of vectors that describe the communication classes of mc.

Given $K$ and $F$, compute the value of deterministic entropy, which is $\sum_t \beta^t x_t' K'K x_t$ with $x_{t+1} = (A - BF + CK) x_t$.

Arguments

• rlq::RBLQ: Instance of RBLQ type
• F::Matrix{Float64} The policy function, a k x n array
• K::Matrix{Float64} The worst case matrix, a j x n array
• x0::Vector{Float64} : The initial condition for state

Returns

• e::Float64 The deterministic entropy

Repeatedly apply a function to search for a fixed point

Approximates $T^∞ v$, where $T$ is an operator (function) and $v$ is an initial guess for the fixed point. Will terminate either when T^{k+1}(v) - T^k v < err_tol or max_iter iterations has been exceeded.

Provided that $T$ is a contraction mapping or similar, the return value will be an approximation to the fixed point of $T$.

Arguments

• T: A function representing the operator $T$
• v::TV: The initial condition. An object of type $TV$
• ;err_tol(1e-3): Stopping tolerance for iterations
• ;max_iter(50): Maximum number of iterations
• ;verbose(2): Level of feedback (0 for no output, 1 for warnings only, 2 for warning and convergence messages during iteration)
• ;print_skip(10) : if verbose == 2, how many iterations to apply between print messages

Returns

• '::TV': The fixed point of the operator $T$. Has type $TV$

Example

using QuantEcon
T(x, μ) = 4.0 * μ * x * (1.0 - x)
x_star = compute_fixed_point(x->T(x, 0.3), 0.4)  # (4μ - 1)/(4μ)

Compute the $v$-greedy policy

Parameters

• ddp::DiscreteDP : Object that contains the model parameters
• ddpr::DPSolveResult : Object that contains result variables

Returns

• sigma::Vector{Int} : Array containing v-greedy policy rule

Notes

modifies ddpr.sigma and ddpr.Tv in place

Compute the $v$-greedy policy.

Arguments

• v::AbstractVector Value function vector of length n
• ddp::DiscreteDP Object that contains the model parameters

Returns

• sigma:: v-greedy policy vector, of lengthn`

computes log-likelihood of entire observations

Arguments

• kn::Kalman: Kalman specifying the model. Initial value must be the prior for t=1 period observation, i.e. $x_{1|0}$.
• y::AbstractMatrix: n x T matrix of observed data. n is the number of observed variables in one period. Each column is a vector of observations at each period.

Returns

• logL::Real: log-likelihood of all observations

Compute and return the optimal state and control sequence, assuming innovation $N(0,1)$

Arguments

• lq::LQ : instance of LQ type
• x0::ScalarOrArray: initial state
• ts_length::Integer(100) : maximum number of periods for which to return process. If lq instance is finite horizon type, the sequenes are returned only for min(ts_length, lq.capT)

Returns

• x_path::Matrix{Float64} : An n x T+1 matrix, where the t-th column represents $x_t$
• u_path::Matrix{Float64} : A k x T matrix, where the t-th column represents $u_t$
• w_path::Matrix{Float64} : A n x T+1 matrix, where the t-th column represents lq.C*N(0,1)

The $D$ operator, mapping $P$ into

\[ D(P) := P + PC(\theta I - C'PC)^{-1} C'P\]

Arguments

• rlq::RBLQ: Instance of RBLQ type
• P::Matrix{Float64} : size is n x n

Returns

• dP::Matrix{Float64} : The matrix $P$ after applying the $D$ operator

Compute a finite-state Markov chain approximation to a VAR(1) process of the form

\[ y_{t+1} = b + By_{t} + \Psi^{\frac{1}{2}}\epsilon_{t+1}\]

where $\epsilon_{t+1}$ is an vector of independent standard normal innovations of length M

P, X = discrete_var(b, B, Psi, Nm, n_moments, method, n_sigmas)

Arguments

• b::Union{Real, AbstractVector} : constant vector of length M. M=1 corresponds scalar case
• B::Union{Real, AbstractMatrix} : M x M matrix of impact coefficients
• Psi::Union{Real, AbstractMatrix} : M x M variance-covariance matrix of the innovations
  • discrete_var only accepts non-singular variance-covariance matrices, Psi.
• Nm::Integer > 3 : Desired number of discrete points in each dimension

Optional

• n_moments::Integer : Desired number of moments to match. The default is 2.
• method::VAREstimationMethod : Specify the method used to determine the grid points. Accepted inputs are Even(), Quantile(), or Quadrature(). Please see the paper for more details.
• n_sigmas::Real : If the Even() option is specified, n_sigmas is used to determine the number of unconditional standard deviations used to set the endpoints of the grid. The default is sqrt(Nm-1).

Returns

• P : Nm^M x Nm^M probability transition matrix. Each row corresponds to a discrete conditional probability distribution over the state M-tuples in X
• X : M x Nm^M matrix of states. Each column corresponds to an M-tuple of values which correspond to the state associated with each row of P

NOTES

• discrete_var only constructs tensor product grids where each dimension contains the same number of points. For this reason it is recommended that this code not be used for problems of more than about 4 or 5 dimensions due to curse of dimensionality issues.
• Future updates will allow for singular variance-covariance matrices and sparse grid specifications.

Reference

• Farmer, L. E., & Toda, A. A. (2017). "Discretizing nonlinear, non‐Gaussian Markov processes with exact conditional moments," Quantitative Economics, 8(2), 651-683.

Given a function f defined on the interval [x1, x2], subdivide the interval into n equally spaced segments, and search for zero crossings of the function. nroot will be set to the number of bracketing pairs found. If it is positive, the arrays xb1[1..nroot] and xb2[1..nroot] will be filled sequentially with any bracketing pairs that are found.

Arguments

• f::Function: The function you want to bracket
• x1::T: Lower border for search interval
• x2::T: Upper border for search interval
• n::Int(50): The number of sub-intervals to divide [x1, x2] into

Returns

• x1b::Vector{T}: Vector of lower borders of bracketing intervals
• x2b::Vector{T}: Vector of upper borders of bracketing intervals

References

This is zbrack from Numerical Recepies Recepies in C++

Approximate the integral of f, given quadrature nodes and weights

Arguments

• f::Function: A callable function that is to be approximated over the domain spanned by nodes.
• nodes::Array: Quadrature nodes
• weights::Array: Quadrature nodes
• args...(Void): additional positional arguments to pass to f
• ;kwargs...(Void): additional keyword arguments to pass to f

Returns

• out::Float64 : The scalar that approximates integral of f on the hypercube formed by [a, b]

Accepts the simulation of a discrete state Markov chain and estimates the transition probabilities

Let $S = s_1, s_2, \ldots, s_N$ with $s_1 < s_2 < \ldots < s_N$ be the discrete states of a Markov chain. Furthermore, let $P$ be the corresponding stochastic transition matrix.

Given a history of observations, $\{X\}_{t=0}^{T}$ with $x_t \in S \forall t$, we would like to estimate the transition probabilities in $P$ with $p_{ij}$ as the ith row and jth column of $P$. For $x_t = s_i$ and $x_{t-1} = s_j$, let $P(x_t | x_{t-1})$ be defined as $p_{i,j}$ element of the stochastic matrix. The likelihood function is then given by

\[ L(\{X\}^t; P) = \text{Prob}(x_1) \prod_{t=2}^{T} P(x_t | x_{t-1})\]

The maximum likelihood estimate is then just given by the number of times a transition from $s_i$ to $s_j$ is observed divided by the number of times $s_i$ was observed.

Note: Because of the estimation procedure used, only states that are observed in the history appear in the estimated Markov chain... It can't divine whether there are unobserved states in the original Markov chain.

For more info, refer to:

• http://www.stat.cmu.edu/~cshalizi/462/lectures/06/markov-mle.pdf
• https://stats.stackexchange.com/questions/47685/calculating-log-likelihood-for-given-mle-markov-chains

Arguments

• X::Vector{T} : Simulated history of Markov states

Returns

• mc::MarkovChain{T} : A Markov chain holding the state values and transition matrix

Given a fixed policy $F$, with the interpretation $u = -F x$, this function computes the matrix $P_F$ and constant $d_F$ associated with discounted cost $J_F(x) = x' P_F x + d_F$.

Arguments

• rlq::RBLQ: Instance of RBLQ type
• F::Matrix{Float64} : The policy function, a k x n array

Returns

• P_F::Matrix{Float64} : Matrix for discounted cost
• d_F::Float64 : Constant for discounted cost
• K_F::Matrix{Float64} : Worst case policy
• O_F::Matrix{Float64} : Matrix for discounted entropy
• o_F::Float64 : Constant for discounted entropy

Method of evaluate_policy that extracts sigma from a DPSolveResult

See other docstring for details

Compute the value of a policy.

Parameters

• ddp::DiscreteDP : Object that contains the model parameters
• sigma::AbstractVector{T<:Integer} : Policy rule vector

Returns

• v_sigma::Array{Float64} : Value vector of sigma, of length n.

Given a function f and an initial guessed range x1 to x2, the routine expands the range geometrically until a root is bracketed by the returned values x1 and x2 (in which case zbrac returns true) or until the range becomes unacceptably large (in which case a ConvergenceError is thrown).

Arguments

• f::Function: The function you want to bracket
• x1::T: Initial guess for lower border of bracket
• x2::T: Initial guess ofr upper border of bracket
• ;ntry::Int(50): The maximum number of expansion iterations
• ;fac::Float64(1.6): Expansion factor (higher ⟶ larger interval size jumps)

Returns

• x1::T: The lower end of an actual bracketing interval
• x2::T: The upper end of an actual bracketing interval

References

This method is zbrac from numerical recipies in C++

Exceptions

• Throws a ConvergenceError if the maximum number of iterations is exceeded

Updates the moments of the time $t$ filtering distribution to the moments of the predictive distribution, which becomes the time $t+1$ prior

Arguments

• k::Kalman An instance of the Kalman filter

Applies Golden-section search to search for the maximum of a function in the interval (a, b)

https://en.wikipedia.org/wiki/Golden-section_search

gridmake(arrays::Union{AbstractVector,AbstractMatrix}...)

Expand one or more vectors (or matrices) into a matrix where rows span the cartesian product of combinations of the input arrays. Each column of the input arrays will correspond to one column of the output matrix. The first array varies the fastest (see example)

Example

julia> x = [1, 2, 3]; y = [10, 20]; z = [100, 200];

julia> gridmake(x, y, z)
12×3 Matrix{Int64}:
 1  10  100
 2  10  100
 3  10  100
 1  20  100
 2  20  100
 3  20  100
 1  10  200
 2  10  200
 3  10  200
 1  20  200
 2  20  200
 3  20  200

gridmake!(out::AbstractMatrix, arrays::AbstractVector...)

Like gridmake, but fills a pre-populated array. out must have size prod(map(length, arrays), dims = length(arrays))

This routine computes the stationary distribution of an irreducible Markov transition matrix (stochastic matrix) or transition rate matrix (generator matrix) $A$.

More generally, given a Metzler matrix (square matrix whose off-diagonal entries are all nonnegative) $A$, this routine solves for a nonzero solution $x$ to $x (A - D) = 0$, where $D$ is the diagonal matrix for which the rows of $A - D$ sum to zero (i.e., $D_{ii} = \sum_j A_{ij}$ for all $i$). One (and only one, up to normalization) nonzero solution exists corresponding to each reccurent class of $A$, and in particular, if $A$ is irreducible, there is a unique solution; when there are more than one solution, the routine returns the solution that contains in its support the first index $i$ such that no path connects $i$ to any index larger than $i$. The solution is normalized so that its 1-norm equals one. This routine implements the Grassmann-Taksar-Heyman (GTH) algorithm (Grassmann, Taksar, and Heyman 1985), a numerically stable variant of Gaussian elimination, where only the off-diagonal entries of $A$ are used as the input data. For a nice exposition of the algorithm, see Stewart (2009), Chapter 10.

Arguments

• A::Matrix{T} : Stochastic matrix or generator matrix. Must be of shape n x n.

Returns

• x::Vector{T} : Stationary distribution of $A$.

References

• W. K. Grassmann, M. I. Taksar and D. P. Heyman, "Regenerative Analysis and Steady State Distributions for Markov Chains, " Operations Research (1985), 1107-1116.
• W. J. Stewart, Probability, Markov Chains, Queues, and Simulation, Princeton University Press, 2009.

This function applies "Hamilton filter" to AbstractVector.

http://econweb.ucsd.edu/~jhamilto/hp.pdf

Arguments

• y::AbstractVector : data to be filtered
• h::Integer : Time horizon that we are likely to predict incorrectly. Original paper recommends 2 for annual data, 8 for quarterly data, 24 for monthly data.
• p::Integer : Number of lags in regression. Must be greater than h.

Note: For seasonal data, it's desirable for p and h to be integer multiples of the number of obsevations in a year. e.g. For quarterly data, h = 8 and p = 4 are recommended.

Returns

• y_cycle::Vector : cyclical component
• y_trend::Vector : trend component

This function applies "Hamilton filter" to <:AbstractVector under random walk assumption.

http://econweb.ucsd.edu/~jhamilto/hp.pdf

Arguments

• y::AbstractVector : data to be filtered
• h::Integer : Time horizon that we are likely to predict incorrectly. Original paper recommends 2 for annual data, 8 for quarterly data, 24 for monthly data.

Note: For seasonal data, it's desirable for h to be an integer multiple of the number of obsevations in a year. e.g. For quarterly data, h = 8 is recommended.

Returns

• y_cycle::Vector : cyclical component
• y_trend::Vector : trend component

apply Hodrick-Prescott filter to AbstractVector.

Arguments

• y::AbstractVector : data to be detrended
• λ::Real : penalty on variation in trend

Returns

• y_cyclical::Vector: cyclical component
• y_trend::Vector: trend component

Get the impulse response corresponding to our model.

Arguments

• arma::ARMA: Instance of ARMA type
• ;impulse_length::Integer(30): Length of horizon for calcluating impulse reponse. Must be at least as long as the p fields of arma

Returns

• psi::Vector{Float64}: psi[j] is the response at lag j of the impulse response. We take psi[1] as unity.

interp(grid::AbstractVector, function_vals::AbstractVector)

Linear interpolation in one dimension

Examples

breaks = cumsum(0.1 .* rand(20))
vals = 0.1 .* sin.(breaks)
li = interp(breaks, vals)

# Do interpolation by treating `li` as a function you can pass scalars to
li(0.2)

# use broadcasting to evaluate at multiple points
li.([0.1, 0.2, 0.3])

Indicate whether the Markov chain mc is aperiodic.

Arguments

• mc::MarkovChain : MarkovChain instance.

Returns

• ::Bool

Indicate whether the Markov chain mc is irreducible.

Arguments

• mc::MarkovChain : MarkovChain instance.

Returns

• ::Bool

is_stable(A)

General function for testing for stability of matrix $A$. Just checks that eigenvalues are less than 1 in absolute value.

Arguments

• A::Matrix The matrix we want to check

Returns

• stable::Bool Whether or not the matrix is stable

Test for stability of linear state space system. First removes the constant row and column.

Arguments

• lss::LSS The linear state space system

Returns

• stable::Bool Whether or not the system is stable

k_array_rank([T=Int], a)

Given an array a of k distinct positive integers, sorted in ascending order, return its ranking in the lexicographic ordering of the descending sequences of the elements, following Combinatorial number system.

Notes

InexactError exception will be thrown, or an incorrect value will be returned without warning if overflow occurs during the computation. It is the user's responsibility to ensure that the rank of the input array fits within the range of T; a sufficient condition for it is binomial(BigInt(a[end]), BigInt(length(a))) <= typemax(T).

Arguments

• T::Type{<:Integer}: The numeric type of ranking to be returned.
• a::Vector{<:Integer}: Array of length k.

Returns

• idx::T: Ranking of a.

A vectorized function that returns the value of the look ahead estimate at the values in the array y.

Arguments

• l::LAE: Instance of LAE type
• y::Array: Array that becomes the y in l.p(l.x, y)

Returns

• psi_vals::Vector: Density at (x, y)

Computes the quadratic sum

\[ V = \sum_{j=0}^{\infty} A^j B A^{j'}\]

$V$ is computed by solving the corresponding discrete lyapunov equation using the doubling algorithm. See the documentation of solve_discrete_lyapunov for more information.

Arguments

• A::Matrix{Float64} : An n x n matrix as described above. We assume in order for convergence that the eigenvalues of $A$ have moduli bounded by unity
• B::Matrix{Float64} : An n x n matrix as described above. We assume in order for convergence that the eigenvalues of $B$ have moduli bounded by unity
• max_it::Int(50) : Maximum number of iterations

Returns

• gamma1::Matrix{Float64} : Represents the value $V$

Create an iterator to calculate the population mean and variance-convariance matrix for both $x_t$ and $y_t$, starting at the initial condition (self.mu_0, self.Sigma_0). Each iteration produces a 4-tuple of items (mu_x, mu_y, Sigma_x, Sigma_y) for the next period.

Arguments

• lss::LSS An instance of the Gaussian linear state space model

Number of states in the Markov chain mc

next_k_array!(a)

Given an array a of k distinct positive integers, sorted in ascending order, return the next k-array in the lexicographic ordering of the descending sequences of the elements, following Combinatorial number system. a is modified in place.

Arguments

• a::Vector{<:Integer}: Array of length k.

Returns

• a::Vector{<:Integer}: View of a.

Examples

julia> n, k = 4, 2;

julia> a = collect(1:2);

julia> while a[end] <= n
           @show a
           next_k_array!(a)
       end
a = [1, 2]
a = [1, 3]
a = [2, 3]
a = [1, 4]
a = [2, 4]
a = [3, 4]

Compute the limit of a Nash linear quadratic dynamic game.

Player i minimizes

\[ \sum_{t=1}^{\infty}(x_t' r_i x_t + 2 x_t' w_i u_{it} +u_{it}' q_i u_{it} + u_{jt}' s_i u_{jt} + 2 u_{jt}' m_i u_{it})\]

subject to the law of motion

\[ x_{t+1} = A x_t + b_1 u_{1t} + b_2 u_{2t}\]

and a perceived control law $u_j(t) = - f_j x_t$ for the other player.

The solution computed in this routine is the $f_i$ and $p_i$ of the associated double optimal linear regulator problem.

Arguments

• A : Corresponds to the above equation, should be of size (n, n)
• B1 : As above, size (n, k_1)
• B2 : As above, size (n, k_2)
• R1 : As above, size (n, n)
• R2 : As above, size (n, n)
• Q1 : As above, size (k_1, k_1)
• Q2 : As above, size (k_2, k_2)
• S1 : As above, size (k_1, k_1)
• S2 : As above, size (k_2, k_2)
• W1 : As above, size (n, k_1)
• W2 : As above, size (n, k_2)
• M1 : As above, size (k_2, k_1)
• M2 : As above, size (k_1, k_2)
• ;beta::Float64(1.0) Discount rate
• ;tol::Float64(1e-8) : Tolerance level for convergence
• ;max_iter::Int(1000) : Maximum number of iterations allowed

Returns

• F1::Matrix{Float64}: (k_1, n) matrix representing feedback law for agent 1
• F2::Matrix{Float64}: (k_2, n) matrix representing feedback law for agent 2
• P1::Matrix{Float64}: (n, n) matrix representing the steady-state solution to the associated discrete matrix ticcati equation for agent 1
• P2::Matrix{Float64}: (n, n) matrix representing the steady-state solution to the associated discrete matrix riccati equation for agent 2

num_compositions(m, n)

The total number of m-part compositions of n, which is equal to (n + m - 1) choose (m - 1).

Arguments

• m::Int : Number of parts of composition
• n::Int : Integer to decompose

Returns

• ::Int : Total number of m-part compositions of n

Updates the moments (cur_x_hat, cur_sigma) of the time $t$ prior to the time $t$ filtering distribution, using current measurement $y_t$. The updates are according to

\[ \hat{x}^F = \hat{x} + \Sigma G' (G \Sigma G' + R)^{-1} (y - G \hat{x}) \\ \Sigma^F = \Sigma - \Sigma G' (G \Sigma G' + R)^{-1} G \Sigma\]

Arguments

• k::Kalman An instance of the Kalman filter
• y The current measurement

Computes nodes and weights for beta distribution.

Arguments

• n::Union{Int, Vector{Int}} : Number of desired nodes along each dimension
• a::Union{Real, Vector{Real}} : First parameter of the beta distribution, along each dimension
• b::Union{Real, Vector{Real}} : Second parameter of the beta distribution, along each dimension

Returns

• nodes::Array{Float64} : An array of quadrature nodes
• weights::Array{Float64} : An array of corresponding quadrature weights

Notes

If any of the parameters to this function are scalars while others are vectors of length n, the the scalar parameter is repeated n times.

References

Miranda, Mario J, and Paul L Fackler. Applied Computational Economics and Finance, MIT Press, 2002.

Computes multivariate Guass-Checbychev quadrature nodes and weights.

Arguments

• n::Union{Int, Vector{Int}} : Number of desired nodes along each dimension
• a::Union{Real, Vector{Real}} : Lower endpoint along each dimension
• b::Union{Real, Vector{Real}} : Upper endpoint along each dimension

Returns

• nodes::Array{Float64} : An array of quadrature nodes
• weights::Array{Float64} : An array of corresponding quadrature weights

Notes

If any of the parameters to this function are scalars while others are vectors of length n, the the scalar parameter is repeated n times.

References

Miranda, Mario J, and Paul L Fackler. Applied Computational Economics and Finance, MIT Press, 2002.

qnwdist(
    d::Distributions.ContinuousUnivariateDistribution, N::Int,
    q0::Real=0.001, qN::Real=0.999, method::Union{T,Type{T}}=Quantile
) where T

Construct N quadrature weights and nodes for distribution d from the quantile q0 to the quantile qN. method can be one of:

• Even: nodes will be evenly spaced between the quantiles
• Quantile: nodes will be placed at evenly spaced quantile values

To construct the weights, consider splitting the nodes into cells centered at each node. Specifically, let notation z_i mean the ith node and let z_{i-1/2} be 1/2 between nodes z_{i-1} and z_i. Then, weights are determined as follows:

• weights[1] = cdf(d, z_{1+1/2})
• weights[N] = 1 - cdf(d, z_{N-1/2})
• weights[i] = cdf(d, z_{i+1/2}) - cdf(d, z_{i-1/2}) for all i in 2:N-1

In effect, this strategy assigns node i all the probability associated with a random variable occuring within the node is cell.

The weights always sum to 1, so they can be used as a proper probability distribution. This means that E[f(x) | x ~ d] ≈ dot(f.(nodes), weights).

Generates equidistributed sequences with property that averages value of integrable function evaluated over the sequence converges to the integral as n goes to infinity.

Arguments

• n::Union{Int, Vector{Int}} : Number of desired nodes along each dimension
• a::Union{Real, Vector{Real}} : Lower endpoint along each dimension
• b::Union{Real, Vector{Real}} : Upper endpoint along each dimension
• kind::AbstractString("N"): One of the following:
  • N - Neiderreiter (default)
  • W - Weyl
  • H - Haber
  • R - pseudo Random

Returns

• nodes::Array{Float64} : An array of quadrature nodes
• weights::Array{Float64} : An array of corresponding quadrature weights

Notes

If any of the parameters to this function are scalars while others are vectors of length n, the the scalar parameter is repeated n times.

References

Miranda, Mario J, and Paul L Fackler. Applied Computational Economics and Finance, MIT Press, 2002.

Computes nodes and weights for beta distribution

Arguments

• n::Union{Int, Vector{Int}} : Number of desired nodes along each dimension
• a::Union{Real, Vector{Real}} : Shape parameter of the gamma distribution, along each dimension. Must be positive. Default is 1
• b::Union{Real, Vector{Real}} : Scale parameter of the gamma distribution, along each dimension. Must be positive. Default is 1

Returns

• nodes::Array{Float64} : An array of quadrature nodes
• weights::Array{Float64} : An array of corresponding quadrature weights

Notes

If any of the parameters to this function are scalars while others are vectors of length n, the the scalar parameter is repeated n times.

References

Miranda, Mario J, and Paul L Fackler. Applied Computational Economics and Finance, MIT Press, 2002.

Computes multivariate Guass-Legendre quadrature nodes and weights.

Arguments

• n::Union{Int, Vector{Int}} : Number of desired nodes along each dimension
• a::Union{Real, Vector{Real}} : Lower endpoint along each dimension
• b::Union{Real, Vector{Real}} : Upper endpoint along each dimension

Returns

• nodes::Array{Float64} : An array of quadrature nodes
• weights::Array{Float64} : An array of corresponding quadrature weights

Notes

If any of the parameters to this function are scalars while others are vectors of length n, the the scalar parameter is repeated n times.

References

Miranda, Mario J, and Paul L Fackler. Applied Computational Economics and Finance, MIT Press, 2002.

Computes quadrature nodes and weights for multivariate uniform distribution

Arguments

• n::Union{Int, Vector{Int}} : Number of desired nodes along each dimension
• mu::Union{Real, Vector{Real}} : Mean along each dimension
• sig2::Union{Real, Vector{Real}, Matrix{Real}}(eye(length(n))) : Covariance structure

Returns

• nodes::Array{Float64} : An array of quadrature nodes
• weights::Array{Float64} : An array of corresponding quadrature weights

Notes

See also the documentation for qnwnorm

References

Miranda, Mario J, and Paul L Fackler. Applied Computational Economics and Finance, MIT Press, 2002.

Computes nodes and weights for multivariate normal distribution.

Arguments

• n::Union{Int, Vector{Int}} : Number of desired nodes along each dimension
• mu::Union{Real, Vector{Real}} : Mean along each dimension
• sig2::Union{Real, Vector{Real}, Matrix{Real}}(eye(length(n))) : Covariance structure

Returns

• nodes::Array{Float64} : An array of quadrature nodes
• weights::Array{Float64} : An array of corresponding quadrature weights

Notes

This function has many methods. I try to describe them here.

n or mu can be a vector or a scalar. If just one is a scalar the other is repeated to match the length of the other. If both are scalars, then the number of repeats is inferred from sig2.

sig2 can be a matrix, vector or scalar. If it is a matrix, it is treated as the covariance matrix. If it is a vector, it is considered the diagonal of a diagonal covariance matrix. If it is a scalar it is repeated along the diagonal as many times as necessary, where the number of repeats is determined by the length of either n and/or mu (which ever is a vector).

If all 3 are scalars, then 1d nodes are computed. mu and sig2 are treated as the mean and variance of a 1d normal distribution

References

Miranda, Mario J, and Paul L Fackler. Applied Computational Economics and Finance, MIT Press, 2002.

Computes multivariate Simpson quadrature nodes and weights.

Arguments

• n::Union{Int, Vector{Int}} : Number of desired nodes along each dimension
• a::Union{Real, Vector{Real}} : Lower endpoint along each dimension
• b::Union{Real, Vector{Real}} : Upper endpoint along each dimension

Returns

• nodes::Array{Float64} : An array of quadrature nodes
• weights::Array{Float64} : An array of corresponding quadrature weights

Notes

If any of the parameters to this function are scalars while others are vectors of length n, the the scalar parameter is repeated n times.

References

Miranda, Mario J, and Paul L Fackler. Applied Computational Economics and Finance, MIT Press, 2002.

Computes multivariate trapezoid quadrature nodes and weights.

Arguments

• n::Union{Int, Vector{Int}} : Number of desired nodes along each dimension
• a::Union{Real, Vector{Real}} : Lower endpoint along each dimension
• b::Union{Real, Vector{Real}} : Upper endpoint along each dimension

Returns

• nodes::Array{Float64} : An array of quadrature nodes
• weights::Array{Float64} : An array of corresponding quadrature weights

Notes

If any of the parameters to this function are scalars while others are vectors of length n, the the scalar parameter is repeated n times.

References

Miranda, Mario J, and Paul L Fackler. Applied Computational Economics and Finance, MIT Press, 2002.

Computes quadrature nodes and weights for multivariate uniform distribution.

Arguments

• n::Union{Int, Vector{Int}} : Number of desired nodes along each dimension
• a::Union{Real, Vector{Real}} : Lower endpoint along each dimension
• b::Union{Real, Vector{Real}} : Upper endpoint along each dimension

Returns

• nodes::Array{Float64} : An array of quadrature nodes
• weights::Array{Float64} : An array of corresponding quadrature weights

Notes

If any of the parameters to this function are scalars while others are vectors of length n, the the scalar parameter is repeated n times.

References

Miranda, Mario J, and Paul L Fackler. Applied Computational Economics and Finance, MIT Press, 2002.

Integrate the d-dimensional function f on a rectangle with lower and upper bound for dimension i defined by a[i] and b[i], respectively; using n[i] points.

Arguments

• f::Function The function to integrate over. This should be a function that accepts as its first argument a matrix representing points along each dimension (each dimension is a column). Other arguments that need to be passed to the function are caught by args... and kwargs...`
• n::Union{Int, Vector{Int}} : Number of desired nodes along each dimension
• a::Union{Real, Vector{Real}} : Lower endpoint along each dimension
• b::Union{Real, Vector{Real}} : Upper endpoint along each dimension
• kind::AbstractString("lege") Specifies which type of integration to perform. Valid values are:
  • "lege" : Gauss-Legendre
  • "cheb" : Gauss-Chebyshev
  • "trap" : trapezoid rule
  • "simp" : Simpson rule
  • "N" : Neiderreiter equidistributed sequence
  • "W" : Weyl equidistributed sequence
  • "H" : Haber equidistributed sequence
  • "R" : Monte Carlo
  • args...(Void): additional positional arguments to pass to f
  • ;kwargs...(Void): additional keyword arguments to pass to f

Returns

• out::Float64 : The scalar that approximates integral of f on the hypercube formed by [a, b]

References

Miranda, Mario J, and Paul L Fackler. Applied Computational Economics and Finance, MIT Press, 2002.

random_discrete_dp([rng], num_states, num_actions[, beta];
                   k=num_states, scale=1)

Generate a DiscreteDP randomly. The reward values are drawn from the normal distribution with mean 0 and standard deviation scale.

Arguments

• rng::AbstractRNG=GLOBAL_RNG : Random number generator.
• num_states::Integer : Number of states.
• num_actions::Integer : Number of actions.
• beta::Real=rand(rng) : Discount factor. Randomly chosen from [0, 1) if not specified.
• ;k::Integer(num_states) : Number of possible next states for each state-action pair. Equal to num_states if not specified.
• scale::Real(1) : Standard deviation of the normal distribution for the reward values.

Returns

• ddp::DiscreteDP : An instance of DiscreteDP.

random_markov_chain([rng], n[, k])

Return a randomly sampled MarkovChain instance with n states, where each state has k states with positive transition probability.

Arguments

• rng::AbstractRNG=GLOBAL_RNG : Random number generator.
• n::Integer : Number of states.
• k::Integer=n : Number of nonzero entries in each column of the matrix. Set to n if none specified.

Returns

• mc::MarkovChain : MarkovChain instance.

Examples

julia> using QuantEcon, Random

julia> rng = MersenneTwister(1234);

julia> mc = random_markov_chain(rng, 3);

julia> mc.p
3×3 LinearAlgebra.Transpose{Float64,Array{Float64,2}}:
 0.590845  0.175952   0.233203
 0.460085  0.106152   0.433763
 0.794026  0.0601209  0.145853

julia> mc = random_markov_chain(rng, 3, 2);

julia> mc.p
3×3 LinearAlgebra.Transpose{Float64,Array{Float64,2}}:
 0.0       0.200586  0.799414
 0.701386  0.0       0.298614
 0.753163  0.246837  0.0

random_stochastic_matrix([rng], n[, k])

Return a randomly sampled n x n stochastic matrix with k nonzero entries for each row.

Arguments

• rng::AbstractRNG=GLOBAL_RNG : Random number generator.
• n::Integer : Number of states.
• k::Integer=n : Number of nonzero entries in each column of the matrix. Set to n if none specified.

Returns

• p::Array : Stochastic matrix.

Find the recurrent classes of the Markov chain mc.

Arguments

• mc::MarkovChain : MarkovChain instance.

Returns

• ::Vector{Vector{Int}} : Vector of vectors that describe the recurrent classes of mc.

Finds the row and column, if any, that correspond to the constant term in a LSS system and removes them to get the matrix that needs to be checked for stability.

Arguments

• lss::LSS The linear state space system

Returns

• A::Matrix The matrix A with constant row and column removed

Simulate num_reps observations of $x_T$ and $y_T$ given $x_0 \sim N(\mu_0, \Sigma_0)$.

Arguments

• lss::LSS An instance of the Gaussian linear state space model.
• t::Int = 10 The period that we want to replicate values for.
• num_reps::Int = 100 The number of replications we want

Returns

• x::Matrix An n x num_reps matrix, where the j-th column is the jth observation of ``xT``
• y::Matrix An k x num_reps matrix, where the j-th column is the jth observation of ``yT``

Find a root of the f on the bracketing inverval [x1, x2] via ridder algo

Arguments

• f::Function: The function you want to bracket
• x1::T: Lower border for search interval
• x2::T: Upper border for search interval
• ;maxiter::Int(500): Maximum number of bisection iterations
• ;xtol::Float64(1e-12): The routine converges when a root is known to lie within xtol of the value return. Should be >= 0. The routine modifies this to take into account the relative precision of doubles.
• ;rtol::Float64(2*eps()):The routine converges when a root is known to lie within rtol times the value returned of the value returned. Should be ≥ 0

Returns

• x::T: The found root

Exceptions

• Throws an ArgumentError if [x1, x2] does not form a bracketing interval
• Throws a ConvergenceError if the maximum number of iterations is exceeded

References

Matches ridder function from scipy/scipy/optimize/Zeros/ridder.c

Solves the robust control problem.

The algorithm here tricks the problem into a stacked LQ problem, as described in chapter 2 of Hansen- Sargent's text "Robustness". The optimal control with observed state is

\[ u_t = - F x_t\]

And the value function is $-x'Px$

Arguments

• rlq::RBLQ: Instance of RBLQ type

Returns

• F::Matrix{Float64} : The optimal control matrix from above
• P::Matrix{Float64} : The positive semi-definite matrix defining the value function
• K::Matrix{Float64} : the worst-case shock matrix $K$, where $w_{t+1} = K x_t$ is the worst case shock

Solve the robust LQ problem

A simple algorithm for computing the robust policy $F$ and the corresponding value function $P$, based around straightforward iteration with the robust Bellman operator. This function is easier to understand but one or two orders of magnitude slower than self.robust_rule(). For more information see the docstring of that method.

Arguments

• rlq::RBLQ: Instance of RBLQ type
• P_init::Matrix{Float64}(zeros(rlq.n, rlq.n)) : The initial guess for the

value function matrix

• ;max_iter::Int(80): Maximum number of iterations that are allowed
• ;tol::Real(1e-8) The tolerance for convergence

Returns

• F::Matrix{Float64} : The optimal control matrix from above
• P::Matrix{Float64} : The positive semi-definite matrix defining the value function
• K::Matrix{Float64} : the worst-case shock matrix $K$, where $w_{t+1} = K x_t$ is the worst case shock

Rouwenhorst's method to approximate AR(1) processes.

The process follows

\[ y_t = \mu + \rho y_{t-1} + \epsilon_t\]

where $\epsilon_t \sim N (0, \sigma^2)$

Arguments

• N::Integer : Number of points in markov process
• ρ::Real : Persistence parameter in AR(1) process
• σ::Real : Standard deviation of random component of AR(1) process
• μ::Real(0.0) : Mean of AR(1) process

Returns

• mc::MarkovChain{Float64} : Markov chain holding the state values and transition matrix

simplex_grid(m, n)

Construct an array consisting of the integer points in the (m-1)-dimensional simplex $\{x \mid x_1 + \cdots + x_m = n, x_i \geq 0\}$, or equivalently, the m-part compositions of n, which are listed in lexicographic order. The total number of the points (hence the length of the output array) is L = (n+m-1)!/(n!*(m-1)!) (i.e., (n+m-1) choose (m-1)).

Arguments

• m::Int : Dimension of each point. Must be a positive integer.
• n::Int : Number which the coordinates of each point sum to. Must be a nonnegative integer.

Returns

• out::Matrix{Int} : Array of shape (m, L) containing the integer points in the simplex, aligned in lexicographic order.

Notes

A grid of the (m-1)-dimensional unit simplex with n subdivisions along each dimension can be obtained by simplex_grid(m, n) / n.

Examples

julia> simplex_grid(3, 4)
3×15 Matrix{Int64}:
 0  0  0  0  0  1  1  1  1  2  2  2  3  3  4
 0  1  2  3  4  0  1  2  3  0  1  2  0  1  0
 4  3  2  1  0  3  2  1  0  2  1  0  1  0  0

References

A. Nijenhuis and H. S. Wilf, Combinatorial Algorithms, Chapter 5, Academic Press, 1978.

simplex_index(x, m, n)

Return the index of the point x in the lexicographic order of the integer points of the (m-1)-dimensional simplex $\{x \mid x_0 + \cdots + x_{m-1} = n\}$.

Arguments

• x::Vector{Int} : Integer point in the simplex, i.e., an array of m nonnegative integers that sum to n.
• m::Int : Dimension of each point. Must be a positive integer.
• n::Int : Number which the coordinates of each point sum to. Must be a nonnegative integer.

Returns

• idx::Int : Index of x.

Fill X with sample paths of the Markov chain mc as columns. The resulting matrix has the state values of mc as elements.

Arguments

• X::Matrix : Preallocated matrix to be filled with sample paths

of the Markov chain mc. The element types in X should be the same as the type of the state values of mc

• mc::MarkovChain : MarkovChain instance.
• ;init=rand(1:n_states(mc)) : Can be one of the following
  • blank: random initial condition for each chain
  • scalar: same initial condition for each chain
  • vector: cycle through the elements, applying each as an initial condition until all columns have an initial condition (allows for more columns than initial conditions)

Simulate one sample path of the Markov chain mc. The resulting vector has the state values of mc as elements.

Arguments

• mc::MarkovChain : MarkovChain instance.
• ts_length::Int : Length of simulation
• ;init::Int=rand(1:n_states(mc)) : Initial state

Returns

• X::Vector : Vector containing the sample path, with length ts_length

Fill X with sample paths of the Markov chain mc as columns. The resulting matrix has the indices of the state values of mc as elements.

Arguments

• X::Matrix{Int} : Preallocated matrix to be filled with indices

of the sample paths of the Markov chain mc.

• mc::MarkovChain : MarkovChain instance.
• ;init=rand(1:n_states(mc)) : Can be one of the following
  • blank: random initial condition for each chain
  • scalar: same initial condition for each chain
  • vector: cycle through the elements, applying each as an initial condition until all columns have an initial condition (allows for more columns than initial conditions)

Simulate one sample path of the Markov chain mc. The resulting vector has the indices of the state values of mc as elements.

Arguments

• mc::MarkovChain : MarkovChain instance.
• ts_length::Int : Length of simulation
• ;init::Int=rand(1:n_states(mc)) : Initial state

Returns

• X::Vector{Int} : Vector containing the sample path, with length ts_length

Compute a simulated sample path assuming Gaussian shocks.

Arguments

• arma::ARMA: Instance of ARMA type
• ;ts_length::Integer(90): Length of simulation
• ;impulse_length::Integer(30): Horizon for calculating impulse response (see also docstring for impulse_response)

Returns

• X::Vector{Float64}: Simulation of the ARMA model arma

Smooth the data in x using convolution with a window of requested size and type.

Arguments

• x::Array: An array containing the data to smooth
• window_len::Int(7): An odd integer giving the length of the window
• window::AbstractString("hanning"): A string giving the window type. Possible values are flat, hanning, hamming, bartlett, or blackman

Returns

• out::Array: The array of smoothed data

Version of smooth where window_len and window are keyword arguments

Arguments

• kn::Kalman: Kalman specifying the model. Initial value must be the prior for t=1 period observation, i.e. $x_{1|0}$.
• y::AbstractMatrix: n x T matrix of observed data. n is the number of observed variables in one period. Each column is a vector of observations at each period.

Returns

• x_smoothed::AbstractMatrix: k x T matrix of smoothed mean of states. k is the number of states.
• logL::Real: log-likelihood of all observations
• sigma_smoothed::AbstractArray k x k x T array of smoothed covariance matrix of states.

Solve the dynamic programming problem.

Parameters

• ddp::DiscreteDP : Object that contains the Model Parameters
• method::Type{T<Algo}(VFI): Type name specifying solution method. Acceptable arguments are VFI for value function iteration or PFI for policy function iteration or MPFI for modified policy function iteration
• ;max_iter::Int(250) : Maximum number of iterations
• ;epsilon::Float64(1e-3) : Value for epsilon-optimality. Only used if method is VFI
• ;k::Int(20) : Number of iterations for partial policy evaluation in modified policy iteration (irrelevant for other methods).

Returns

• ddpr::DPSolveResult{Algo} : Optimization result represented as a DPSolveResult. See DPSolveResult for details.

Solves the discrete lyapunov equation.

The problem is given by

\[ AXA' - X + B = 0\]

$X$ is computed by using a doubling algorithm. In particular, we iterate to convergence on $X_j$ with the following recursions for $j = 1, 2, \ldots$ starting from $X_0 = B, a_0 = A$:

\[ a_j = a_{j-1} a_{j-1} \\ X_j = X_{j-1} + a_{j-1} X_{j-1} a_{j-1}'\]

Arguments

• A::Matrix{Float64} : An n x n matrix as described above. We assume in order for convergence that the eigenvalues of $A$ have moduli bounded by unity
• B::Matrix{Float64} : An n x n matrix as described above. We assume in order for convergence that the eigenvalues of $B$ have moduli bounded by unity
• max_it::Int(50) : Maximum number of iterations

Returns

• gamma1::Matrix{Float64} Represents the value $X$

Solves the discrete-time algebraic Riccati equation

The prolem is defined as

\[ X = A'XA - (N + B'XA)'(B'XB + R)^{-1}(N + B'XA) + Q\]

via a modified structured doubling algorithm. An explanation of the algorithm can be found in the reference below.

Arguments

• A : k x k array.
• B : k x n array
• R : n x n, should be symmetric and positive definite
• Q : k x k, should be symmetric and non-negative definite
• N::Matrix{Float64}(zeros(size(R, 1), size(Q, 1))) : n x k array
• tolerance::Float64(1e-10) Tolerance level for convergence
• max_iter::Int(50) : The maximum number of iterations allowed

Note that A, B, R, Q can either be real (i.e. k, n = 1) or matrices.

Returns

• X::Matrix{Float64} The fixed point of the Riccati equation; a k x k array representing the approximate solution

References

Chiang, Chun-Yueh, Hung-Yuan Fan, and Wen-Wei Lin. "STRUCTURED DOUBLING ALGORITHM FOR DISCRETE-TIME ALGEBRAIC RICCATI EQUATIONS WITH SINGULAR CONTROL WEIGHTING MATRICES." Taiwanese Journal of Mathematics 14, no. 3A (2010): pp-935.

Compute the spectral density function.

The spectral density is the discrete time Fourier transform of the autocovariance function. In particular,

\[ f(w) = \sum_k \gamma(k) \exp(-ikw)\]

where $\gamma$ is the autocovariance function and the sum is over the set of all integers.

Arguments

• arma::ARMA: Instance of ARMA type
• ;two_pi::Bool(true): Compute the spectral density function over $[0, \pi]$ if false and $[0, 2 \pi]$ otherwise.
• ;res(1200) : If res is a scalar then the spectral density is computed at res frequencies evenly spaced around the unit circle, but if res is an array then the function computes the response at the frequencies given by the array

Returns

• w::Vector{Float64}: The normalized frequencies at which h was computed, in radians/sample
• spect::Vector{Float64} : The frequency response

Compute stationary distributions of the Markov chain mc, one for each recurrent class.

Arguments

• mc::MarkovChain{T} : MarkovChain instance.

Returns

• stationary_dists::Vector{Vector{T1}} : Vector of vectors that represent stationary distributions, where the element type T1 is Rational if T is Int (and equal to T otherwise).

Compute the moments of the stationary distributions of $x_t$ and $y_t$ if possible. Computation is by iteration, starting from the initial conditions lss.mu_0 and lss.Sigma_0

Arguments

• lss::LSS An instance of the Guassian linear state space model
• ;max_iter::Int = 200 The maximum number of iterations allowed
• ;tol::Float64 = 1e-5 The tolerance level one wishes to achieve

Returns

• mu_x::Vector Represents the stationary mean of $x_t$
• mu_y::Vector Represents the stationary mean of $y_t$
• Sigma_x::Matrix Represents the var-cov matrix
• Sigma_y::Matrix Represents the var-cov matrix

Computes value and policy functions in infinite horizon model.

Arguments

• lq::LQ : instance of LQ type

Returns

• P::ScalarOrArray : n x n matrix in value function representation $V(x) = x'Px + d$
• d::Real : Constant in value function representation
• F::ScalarOrArray : Policy rule that specifies optimal control in each period

Notes

This function updates the P, d, and F fields on the lq instance in addition to returning them

Non-mutating routine for solving for P, d, and F in infinite horizon model

See docstring for stationary_values! for more explanation

Tauchen's (1996) method for approximating AR(1) process with finite markov chain

The process follows

\[ y_t = \mu + \rho y_{t-1} + \epsilon_t\]

where $\epsilon_t \sim N (0, \sigma^2)$

Arguments

• N::Integer: Number of points in markov process
• ρ::Real : Persistence parameter in AR(1) process
• σ::Real : Standard deviation of random component of AR(1) process
• μ::Real(0.0) : Mean of AR(1) process
• n_std::Real(3) : The number of standard deviations to each side the process should span

Returns

• mc::MarkovChain : Markov chain holding the state values and transition matrix

Updates cur_x_hat and cur_sigma given array y of length k. The full update, from one period to the next

Arguments

• k::Kalman An instance of the Kalman filter
• y An array representing the current measurement

Update P and d from the value function representation in finite horizon case

Arguments

• lq::LQ : instance of LQ type

Returns

• P::ScalarOrArray : n x n matrix in value function representation $V(x) = x'Px + d$
• d::Real : Constant in value function representation

Notes

This function updates the P and d fields on the lq instance in addition to returning them

Computes the expected discounted quadratic sum

\[ q(x_0) = \mathbb{E} \sum_{t=0}^{\infty} \beta^t x_t' H x_t\]

Here ${x_t}$ is the VAR process $x_{t+1} = A x_t + C w_t$ with ${w_t}$ standard normal and $x_0$ the initial condition.

Arguments

• A::Union{Float64, Matrix{Float64}} The n x n matrix described above (scalar) if n = 1
• C::Union{Float64, Matrix{Float64}} The n x n matrix described above (scalar) if n = 1
• H::Union{Float64, Matrix{Float64}} The n x n matrix described above (scalar) if n = 1
• beta::Float64: Discount factor in (0, 1)
• x_0::Union{Float64, Vector{Float64}} The initial condtion. A conformable array (of length n) or a scalar if n = 1

Returns

• q0::Float64 : Represents the value $q(x_0)$

Notes

The formula for computing $q(x_0)$ is $q(x_0) = x_0' Q x_0 + v$ where

• $Q$ is the solution to $Q = H + \beta A' Q A$ and
• $v = \frac{trace(C' Q C) \beta}{1 - \beta}$

@def_sim sim_name default_type_params begin
    obs_typedef
end

Given a type definition for a single observation in a simulation (obs_typedef), evaluate that type definition as is, but also creates a second type named sim_name as well as various methods on the new type.

The fields of sim_name will have the same name as the fields of obs_typedef, but will be arrays of whatever the type of the corresponding obs_typedef field was. The intention is for sim_name to be a struct of arrays (see https://en.wikipedia.org/wiki/AOSandSOA). If you want an array of structs, just simply collect an array of instances of the type defined in obs_typedef. The struct of arrays storage format has better cache efficiency and data locality if you want to operate on all values of a particular field at once, rather than all the fields of a particular value.

In addition to the new type sim_name, the following methods will be defined:

• sim_name(sz::NTuple{N,Int}). This is a constructor for sim_name that allocates arrays of size sz for each field. If obs_typedef inlcuded any type parameters, then the default values (specified in default_type_params) will be used.
• Base.endof(::sim_name): equal to the length of any of its fields
• Base.length(::sim_name): equal to the length of any of its fields
• The iterator protocol for sim_name. The type of each element of the iterator is the type defined in obs_typedef. This amounts tho defining the following methods
  • Base.start(::sim_name)::Int
  • Base.next(::sim_name, ::Int)::Tuple{Observation,Int}
  • Base.done(::sim_name, ::Int)::Bool
• Base.getindex(sim::sim_name, ix::Int). This implements linear indexing for sim_name and will return an instance of the type defined in obs_typedef

Example

NOTE: the using MacroTools and call to MacroTools.prettify is not necessary and is only used here to clean up the output so it is easier to read

julia> using MacroTools

julia> macroexpand(:(@def_sim Simulation (T => Float64,) struct Observation{T<:Number}
           c::T
           k::T
           i_z::Int
       end
       )) |> MacroTools.prettify
quote
    struct Simulation{prairiedog, T <: Number}
        c::Array{T, prairiedog}
        k::Array{T, prairiedog}
        i_z::Array{Int, prairiedog}
    end
    function Simulation{prairiedog}(sz::NTuple{prairiedog, Int})
        c = Array{Float64, prairiedog}(sz)
        k = Array{Float64, prairiedog}(sz)
        i_z = Array{Int, prairiedog}(sz)
        Simulation(c, k, i_z)
    end
    struct Observation{T <: Number}
        c::T
        k::T
        i_z::Int
    end
    Base.endof(sim::Simulation) = length(sim.c)
    Base.length(sim::Simulation) = endof(sim)
    Base.start(sim::Simulation) = 1
    Base.next(sim::Simulation, ix::Int) = (sim[ix], ix + 1)
    Base.done(sim::Simulation, ix::Int) = ix >= length(sim)
    function Base.getindex(sim::Simulation, ix::Int)
        $(Expr(:boundscheck, true))
        if ix > length(sim)
            throw(BoundsError("$(length(sim))-element Simulation at index $(ix)"))
        end
        $(Expr(:boundscheck, :pop))
        $(Expr(:inbounds, true))
        out = Observation(sim.c[ix], sim.k[ix], sim.i_z[ix])
        $(Expr(:inbounds, :pop))
        return out
    end
end

DPSolveResult is an object for retaining results and associated metadata after solving the model

Parameters

• ddp::DiscreteDP : DiscreteDP object

Returns

• ddpr::DPSolveResult : DiscreteDP results object

types specifying the method for discrete_var

Define Matrix Multiplication between 3-dimensional matrix and a vector

Matrix multiplication over the last dimension of $A$

Make multiple draws from the discrete distribution represented by a DiscreteRV instance

Arguments

• d::DiscreteRV: The DiscreteRV type representing the distribution
• k::Int

Returns

• out::Vector{Int}: k draws from d

Make a single draw from the discrete distribution.

Arguments

• d::DiscreteRV: The DiscreteRV type represetning the distribution

Returns

• out::Int: One draw from the discrete distribution

Private method implementing compute_sequence when state is a scalar

Private method implementing compute_sequence when state is a scalar

Generate a_indptr; stored in out. s_indices is assumed to be in sorted order.

Parameters

• num_states::Integer
• s_indices::AbstractVector{T}
• out::AbstractVector{T} : with length = num_states + 1

Check whether s_indices and a_indices are sorted in lexicographic order.

Parameters

s_indices, a_indices : Vectors

Returns

bool: Whether s_indices and a_indices are sorted.

_random_stochastic_matrix([rng], n, m; k=n)

Generate a "non-square column stochstic matrix" of shape (n, m), which contains as columns m probability vectors of length n with k nonzero entries.

Arguments

• rng::AbstractRNG=GLOBAL_RNG : Random number generator.
• n::Integer : Number of states.
• m::Integer : Number of probability vectors.
• ;k::Integer(n) : Number of nonzero entries in each column of the matrix. Set to n if none specified.

Returns

• p::Array : Array of shape (n, m) containing m probability vectors of length n as columns.

Modified Policy Function Iteration

Policy Function Iteration

NOTE: The epsilon is ignored in this method. It is only here so dispatch can go from solve(::DiscreteDP, ::Type{Algo}) to any of the algorithms. See solve for further details

Impliments Value Iteration NOTE: See solve for further details

Return combinations of each column of matrix A. It is simiplifying allcomb2 by using gridmake from QuantEcon

Arguments

• A::AbstractMatrix : N x M Matrix

Returns

• N^M x M Matrix, combination of each row of A.

Example

julia> allcomb3([1 4 7;
                 2 5 8;
                 3 6 9]) # numerical input
27×3 Array{Int64,2}:
 1  4  7
 1  4  8
 1  4  9
 1  5  7
 1  5  8
 1  5  9
 1  6  7
 1  6  8
 1  6  9
 2  4  7
 ⋮
 2  6  9
 3  4  7
 3  4  8
 3  4  9
 3  5  7
 3  5  8
 3  5  9
 3  6  7
 3  6  8
 3  6  9

construct one-dimensional quantile grid of states

Argument

• Sigma::AbstractMatrix : variance-covariance matrix of the standardized process
• Nm::Integer : number of grid points
• M::Integer : number of variables (M=1 corresponds to AR(1))
• n_sigmas::Real : number of standard error determining end points of grid
• method::Quntile : method for grid making

Return

• y1D : M x Nm matrix of variable grid
• y1Dbounds : bounds of each grid bin

construct one-dimensional quadrature grid of states

Argument

• ::ScalarOrArray : not used
• Nm::Integer : number of grid points
• M::Integer : number of variables (M=1 corresponds to AR(1))
• n_sigmas::Real : not used
• method::Quadrature : method for grid making

Return

• y1D : M x Nm matrix of variable grid
• weights : weights on each grid

construct one-dimensional evenly spaced grid of states

Argument

• Sigma::ScalarOrArray : variance-covariance matrix of the standardized process
• Nm::Integer : number of grid points
• M::Integer : number of variables (M=1 corresponds to AR(1))
• n_sigmas::Real : number of standard error determining end points of grid
• method::Even : method for grid making

Return

• y1D : M x Nm matrix of variable grid
• nothing : nothing of type Void

construct prior guess for quantile grid method

Arguments

• cond_mean::AbstractVector : conditional Mean of each variable
• Nm::Integer : number of grid points
• ::AbstractMatrix : grid of variable
• y1Dbounds::AbstractMatrix : bounds of each grid bin
• method::Quantile : method for grid making

construct prior guess for quadrature grid method

Arguments

• cond_mean::AbstractVector : conditional Mean of each variable
• Nm::Integer : number of grid points
• y1D::AbstractMatrix : grid of variable
• weights::AbstractVector : weights of grid y1D
• method::Quadrature : method for grid making

construct prior guess for evenly spaced grid method

Arguments

• cond_mean::AbstractVector : conditional Mean of each variable
• Nm::Integer : number of grid points
• y1D::AbstractMatrix : grid of variable
• ::AbstractMatrix : bounds of each grid bin
• method::Even : method for grid making

Compute a discrete state approximation to a distribution with known moments, using the maximum entropy procedure proposed in Tanaka and Toda (2013)

p, lambda_bar, moment_error = discrete_approximation(D, T, Tbar, q, lambda0)

Arguments

• D::AbstractVector : vector of grid points of length N. N is the number of points at which an approximation is to be constructed.
• T::Function : A function that accepts a single AbstractVector of length N and returns an L x N matrix of moments evaluated at each grid point, where L is the number of moments to be matched.
• Tbar::AbstractVector : length L vector of moments of the underlying distribution which should be matched

Optional

• q::AbstractVector : length N vector of prior weights for each point in D. The default is for each point to have an equal weight.
• lambda0::AbstractVector : length L vector of initial guesses for the dual problem variables. The default is a vector of zeros.

Returns

• p : (1 x N) vector of probabilties assigned to each grid point in D.
• lambda_bar : length L vector of dual problem variables which solve the maximum entropy problem
• moment_error : vector of errors in moments (defined by moments of discretization minus actual moments) of length L

Compute gradient of objective function

Returns

• grad : length L gradient vector of the objective function evaluated at lambda

Compute hessian of objective function

Returns

• hess : L x L hessian matrix of the objective function evaluated at lambda

Compute the maximum entropy objective function used in discrete_approximation

obj = entropy_obj(lambda, Tx, Tbar, q)

Arguments

• lambda::AbstractVector : length L vector of values of the dual problem variables
• Tx::AbstractMatrix : L x N matrix of moments evaluated at the grid points specified in discrete_approximation
• Tbar::AbstractVector : length L vector of moments of the underlying distribution which should be matched
• q::AbstractVector : length N vector of prior weights for each point in the grid.

Returns

• obj : scalar value of objective function evaluated at lambda

fix(x)

Round x towards zero. For arrays there is a mutating version fix!

Simple method to return an element $Z$ in the Riccati equation solver whose type is Float64 (to be accepted by the cond() function)

Arguments

• BB::Float64 : result of $B' B$
• gamma::Float64 : parameter in the Riccati equation solver
• R::Float64

Returns

• ::Float64 : element $Z$ in the Riccati equation solver

Simple method to return an element $Z$ in the Riccati equation solver whose type is Float64 (to be accepted by the cond() function)

Arguments

• BB::Union{Vector, Matrix} : result of $B' B$
• gamma::Float64 : parameter in the Riccati equation solver
• R::Float64

Returns

• ::Float64 : element $Z$ in the Riccati equation solver

Simple method to return an element $Z$ in the Riccati equation solver whose type is Matrix (to be accepted by the cond() function)

Arguments

• BB::Matrix : result of $B' B$
• gamma::Float64 : parameter in the Riccati equation solver
• R::Matrix

Returns

• ::Matrix : element $Z$ in the Riccati equation solver

Arguments

• kn::Kalman: Kalman specifying the model.
• x_fi::Vector: filtered mean of state for period $t$
• sigma_fi::Matrix: filtered covariance matrix of state for period $t$
• sigma_fo::Matrix: forecast of covariance matrix of state for period $t+1$ conditional on period $t$ observations
• x_s1::Vector: smoothed mean of state for period $t+1$
• sigma_s1::Matrix: smoothed covariance of state for period $t+1$

Returns

• x_s1::Vector: smoothed mean of state for period $t$
• sigma_s1::Matrix: smoothed covariance of state for period $t$

Same as gth_solve, but overwrite the input A, instead of creating a copy.

computes log-likelihood of period $t$

Arguments

• kn::Kalman: Kalman specifying the model. Current values must be the forecast for period $t$ observation conditional on $t-1$ observation.
• y::AbstractVector: Respondentbservations at period $t$

Returns

• logL::Real: log-likelihood of observations at period $t$

find a unitary matrix U such that the diagonal components of U'AU is as close to a multiple of identity matrix as possible

Arguments

• A::AbstractMatrix : square matrix

Returns

• U : unitary matrix
• fval : minimum value

Compute the moment defining function used in discrete_approximation

T = polynomial_moment(X, mu, scaling_factor, mMoments)

Arguments:

• X::AbstractVector : length N vector of grid points
• mu::Real : location parameter (conditional mean)
• scaling_factor::Real : scaling factor for numerical stability. (typically largest grid point)
• n_moments::Integer : number of polynomial moments

Return

• T : moment defining function used in discrete_approximation

random_probvec([rng], k[, m])

Return m randomly sampled probability vectors of size k.

Arguments

• rng::AbstractRNG=GLOBAL_RNG : Random number generator.
• k::Integer : Size of each probability vector.
• m::Integer : Number of probability vectors.

Returns

• a::Array : Matrix of shape (k, m), or Vector of shape (k,) if m is not specified, containing probability vector(s) as column(s).

Populate out with max_a vals(s, a), where vals is represented as a Vector of size (num_sa_pairs,).

Populate out with max_a vals(s, a), where vals is represented as a Vector of size (num_sa_pairs,).

Also fills out_argmax with the cartesiean index associated with the indmax in each row

Populate out with max_a vals(s, a), where vals is represented as a AbstractMatrix of size (num_states, num_actions).

Also fills out_argmax with the column number associated with the indmax in each row

Populate out with max_a vals(s, a), where vals is represented as a AbstractMatrix of size (num_states, num_actions).

Return the Vector max_a vals(s, a), where vals is represented as a AbstractMatrix of size (num_states, num_actions).

return standerdized VAR(1) representation

Arguments

• b::AbstractVector : M x 1 constant term vector
• B::AbstractMatrix : M x M matrix of impact coefficients
• Psi::AbstractMatrix : M x M variance-covariance matrix of innovations
• M::Intger : number of variables of the VAR(1) model

Returns

• A::Matirx : impact coefficients of standardized VAR(1) process
• C::AbstractMatrix : variance-covariance matrix of standardized model innovations
• mu::AbstractVector : mean of the standardized VAR(1) process
• Sigma::AbstractMatrix : variance-covariance matrix of the standardized VAR(1) process

return standerdized AR(1) representation

Arguments

• b::Real : constant term
• B::Real : impact coefficient
• Psi::Real : variance of innovation
• M::Integer == 1 : must be one since the function is for AR(1)

Returns

• A::Real : impact coefficient of standardized AR(1) process
• C::Real : standard deviation of the innovation
• mu::Real : mean of the standardized AR(1) process
• Sigma::Real : variance of the standardized AR(1) process

If A is already dense, return A as is

Custom version of full, which allows convertion to type T

check persistency when method is Quadrature and give warning if needed

Arguments

• B::Union{Real, AbstractMatrix} : impact coefficient
• method::VAREstimationMethod : method for grid making

Returns

• nothing