Uncertain Expansion Example Notebook

Uncertain Expansion Example Notebook #

In Chapter 11, we described the small-noise expansion method and the recursive utility framework. This section covers the usage of Expansion Suite, an open source Python toolbox which uses these methods to solve discrete nonlinear DSGE models. We use an example model taken from Hansen, Khorrami and Tourre (2024).

This notebook provides both written explanations and accompanying code. The simplest way to run the code is to open the Google Colab link at the top of this page. Otherwise, the reader can copy the code snippets onto their local machine.

\(\newcommand{\eqdef}{\stackrel{\text{def}}{=}}\)

1 Model Setup #

We employ a discrete-time approximation of the single-capital model in Section 4.4 of [Hansen et al., 2024], as outlined in Section 11.7.1.

Let \({\widehat K}_t = \log K_t,\), \({\widehat C}_t = \log C_t\), \(\widehat{Z}_t^2 = \log Z_t^2\) and \(D_t = \frac{I_t}{K_t}\). The model uses three states: \(X_t = \begin{bmatrix}{\widehat K}_{t+\epsilon}-{\widehat K}_{t}, Z_t^1, \widehat{Z}_t^2\end{bmatrix}\), as well as two jump variables \(J_t = \begin{bmatrix}{\widehat C}_t, D_t \end{bmatrix}\).

The law of motions for the three states are:

\[\begin{align} {{\widehat K}_{t+\epsilon}} - {\widehat K}_t = \hspace{.2cm} & \epsilon \left[ {\frac 1 \zeta} \log \left( 1 + \zeta \frac{I_t}{K_t} \right) + \nu_k Z_{1,t} - \iota_k \right] \cr \hspace{.2cm} & - \frac{\epsilon}{2}|\sigma_k|^2 \exp\left( {Z}_{2,t} \right) + \sqrt{\epsilon} \exp\left(\frac 1 2 { Z}_{2,t} \right) {\sigma}_k W_{t+\epsilon} \end{align} \]

\[ \begin{align} {Z}_{2,t+\epsilon} - {Z}_{2, t} = \hspace{.2cm} & - \epsilon \nu_2 \left[ 1 - {\mu_2} \exp\left( - {Z}_{2,t} \right) \right] \cr \hspace{.2cm}& - {\frac \epsilon 2} |\sigma_2|^2 \exp\left( - {Z}_{2,t} \right) + \sqrt{\epsilon} \exp\left( - {\frac 1 2} {Z}_{2,t} \right) \sigma_2 W_{t+\epsilon}. \end{align} \]

We can solve for the jump variables using the first-order condition for \(D_t\) and the output constraint:

\[ {\widehat C}_t = \log \left( \alpha - D_t \right) + {\widehat K}_t \]

\[\begin{split} \begin{align*} 0 &= - (1-\beta) \exp\left[(1- \rho) \left({\widehat C}_t - {\widehat G}_t \right)\right] \left( \frac{1}{\alpha - D_t} \right) \\&+ \beta \exp\left[(1 - \rho) \left({\widehat R}_t - {\widehat G}_t \right)\right] \left(\frac{1}{1 + \zeta D_t} \right) \end{align*} \end{split}\]

1.2 Rewriting First-Order Conditions #

To input this model into our code, we need to write it in the form introduced in Section 11.9, \(Q_t {\mathbb E} \left( N_{t+\epsilon} H_{t+\epsilon} \mid {\mathfrak A}_t \right) + P_tL_{t} - M_t = 0\), where

\[ \begin{align*} Q_t \eqdef & \hspace{.2cm} \beta \exp\left[(1-\rho) \left( {\widehat R}_t - {\widehat V}_t \right) \right] \cr \cr P_t \eqdef & \hspace{.2cm} (1-\beta) \exp\left[(\rho-1) \left( {\widehat V}_t - {\widehat G}_t \right) \right] \cr \cr H_{t+\epsilon} \eqdef & \hspace{.2cm} \begin{bmatrix} \psi_{d'}^x (D_t, X_t, W_{t+\epsilon})' & \psi_{d'}^g (D_t, X_t, W_{t+\epsilon}) \cr \psi_{x'}^x (D_t, X_t, W_{t+\epsilon})' & \psi_{g}^x (D_t, X_t, W_{t+\epsilon})' \cr \psi_{g'}^x (D_t, X_t, W_{t+\epsilon})' & 1 \end{bmatrix} \begin{bmatrix} MX_{t+\epsilon} \cr MG_{t+\epsilon} \end{bmatrix} \ \cr \cr L_t \eqdef & \hspace{.2cm} \begin{bmatrix} \exp \left[ (1-\rho) \left({\widehat C}_t - {\widehat G}_t \right) \right] \begin{bmatrix} \kappa_d(D_t,X_t) \cr \kappa_x(D_t, X_t) \cr 1 \end{bmatrix} \end{bmatrix} \cr M_t \eqdef & \hspace{.2cm} \begin{bmatrix} 0 \cr MG_t \cr MX_t \end{bmatrix} . \end{align*} \]

By noting that \(MG_t=1\), we can express the model above simply as:

\[\begin{split} H_{t+\epsilon} = \begin{bmatrix} 0 \\ \frac \epsilon {1+\zeta D_t} \end{bmatrix} \end{split}\]

\[\begin{split} L_{t} = \begin{bmatrix} 0 \\ \exp\left[(1-\rho)(\widehat{C}_t-\widehat{G}_t)\right] \frac 1 {\alpha-D_t} \end{bmatrix} \end{split}\]

\[\begin{split} U_t = \begin{bmatrix} \exp{(\widehat{C}_t - \widehat{K}_t)} + D_t - \alpha \\ 0 \end{bmatrix} \end{split}\]

We stack the system above with the state law of motions and solve the system given the explicit approximations for \(Q_t\) and \(N_{t+\epsilon}\) detailed in the appendix of Chapter 11.

2 Inputs #

The Expansion Suite uses the function uncertain_expansion to approximate a solution to the above system locally. It takes the following inputs:

The output is of class ModelSolution, which stores the coefficients for the linear-quadratic approximation for the jump and state variables; continuation values; consumption growth; and log change of measure, as well as the steady-state values of each variables.

We will now walk through the computation using the example above.

2.1 Libraries #

Begin by installing the following libraries and downloading RiskUncertaintyValue, which contains the functions required to solve the model:

import os
import sys
workdir = os.getcwd()
!git clone https://github.com/lphansen/RiskUncertaintyValue 
workdir = os.getcwd() + '/RiskUncertaintyValue'             
sys.path.insert(0, workdir+'/src')                        
import numpy as np
import autograd.numpy as anp
from scipy import optimize
np.set_printoptions(suppress=True)
np.set_printoptions(linewidth=200)
from IPython.display import display, HTML
from BY_example_sol import disp
display(HTML("<style>.container { width:97% !important; }</style>"))
import warnings
warnings.filterwarnings("ignore")
import matplotlib.pyplot as plt

from lin_quad import LinQuadVar
from lin_quad_util import E, concat, next_period, cal_E_ww, lq_sum, N_tilde_measure, E_N_tp1
from uncertain_expansion import uncertain_expansion
np.set_printoptions(suppress=True)
import pickle
from scipy.optimize import fsolve

fatal: destination path 'RiskUncertaintyValue' already exists and is not an empty directory.

2.2 Equilibrium Condition Function #

The first input for uncertain_expansion requires us to define the equilibrium conditions for the model. Therefore we define a function eq which takes in a list of jump and state variables as inputs, and outputs the equilibrium equations in the form of equations 1 and 2. The inputs required:

Variable	Type	Description	Corresponds to
`Var_t`	array-like	Vector of jump and state variables in the current period	\((J_t, X_t)\)
`Var_tp1`	array-like	Vector of jump and state variables in the next period	\((J_{t+\epsilon}, X_{t+\epsilon})\)
`W_tp1`	array-like	Vector of shock variables in the next period	\((W_{t+\epsilon})\)
`q`	float	Perturbation parameter	\(q\)
`mode`	string	By default, this function returns the LHS of equation (1) and (2). Set mode to ‘psi1’, ‘psi2’, or ‘phi’ to return the corresponding function in the equilibrium conditions
`recursive_ss`	array-like	Steady-state values of \(N\) and \(Q\)
`args`	tuple of float/ndarray	Preference and model parameters

While writing out eq, ensure that:

State variables follow jump variables.
The number of jump variables equals the number of forward-looking equilibrium conditions.
The number of state variables equals to the number of state evolution equations.
All values are expressed as type float. Otherwise, this may cause type errors.

For the example model, we write:

def eq_onecap_3d(Var_t, Var_tp1, W_tp1, q, mode, recursive_ss, *args):
    
    # Unpack model parameters
    γ, β, ρ, a, ϕ_1, ϕ_2, α_k, beta_k, beta_z, sigma_k1, sigma_k2, sigma_k3, sigma_z1, sigma_z2, sigma_z3, sigma_y1, sigma_y2, sigma_y3, beta_2, mu_2, ε = args
    w1_tp1, w2_tp1, w3_tp1 = W_tp1.ravel()

    # Unpack model variables
    log_cmk_t, log_imk_t, gk_t, Z_t, Y_t = Var_t.ravel()
    log_cmk_tp1, log_imk_tp1, gk_tp1, Z_tp1, Y_tp1 = Var_tp1.ravel() 

    imk_t = anp.exp(log_imk_t)
    cmk_t = anp.exp(log_cmk_t)

    q_t = β*anp.exp((1.-ρ)*recursive_ss[0])
    p_t = (1-β)*anp.exp((ρ-1.)*recursive_ss[1])

    # State evolution processes
    dgk = ε * (ϕ_1 * anp.log(1.+ϕ_2*imk_t)-α_k + beta_k*Z_t - q**2* 0.5 * (sigma_k1**2 + sigma_k2**2 + sigma_k3**2) * anp.exp(Y_t)) +  anp.sqrt(ε) *(anp.exp(0.5*Y_t))*(sigma_k1*w1_tp1+sigma_k2*w2_tp1+sigma_k3*w3_tp1)
    # Equation for states at time t
    phi_1 = gk_tp1 - dgk
    phi_2 = Z_tp1 - Z_t + ε * (beta_z*Z_t) - anp.sqrt(ε)*(anp.exp(0.5*Y_t))*(sigma_z1*w1_tp1 + sigma_z2*w2_tp1+ sigma_z3*w3_tp1)
    phi_3 = Y_tp1 - Y_t + ε * beta_2 * (1 - mu_2 * anp.exp(-Y_t)) \
        + q**2 * 0.5 * (sigma_y1**2 + sigma_y2**2 + sigma_y3**2) * anp.exp(-Y_t) * ε \
        - anp.exp(-0.5 * Y_t) * (sigma_y1*w1_tp1+sigma_y2*w2_tp1+sigma_y3*w3_tp1) * anp.sqrt(ε)

    ## Output constraint
    qnh_1 = 0.
    pl_1 = 0.
    u_1 = -a + cmk_t + imk_t

    ## Investment FOC
    qnh_2 =  ε/(1.+ϕ_2*imk_t) 
    pl_2 = - anp.exp((1.-ρ)*log_cmk_t) * 1/(a-imk_t)
    u_2 = 0.

    if mode == 'H':
        return anp.array([qnh_1, qnh_2]) 
    
    if mode == 'L':
        return anp.array([pl_1, pl_2]) 

    return anp.array([
        q_t*qnh_1 - u_1,
        q_t*qnh_2 + p_t* pl_2 - u_2,
        phi_1, phi_2, phi_3])

2.3 Steady State Function #

This function takes model parameters as input and outputs the steady states for each variable following the variable order specified in Var_t and Var_tp1. To find the steady state, we remove the time subscripts from each variable and solve the system of equations.

def ss_onecap_3d(*args, return_recursive=False):
    γ, β, ρ, a, ϕ_1, ϕ_2, α_k, beta_k, beta_z, sigma_k1, sigma_k2, sigma_k3, sigma_z1, sigma_z2, sigma_z3, sigma_y1, sigma_y2, sigma_y3, beta_2, mu_2, ε = args

    def f(x):
        vmk, imk, log_cmk, log_gk = x

        recursive = anp.exp((1 - ρ) * vmk) - (1 - β) * anp.exp((1 - ρ) * log_cmk) - β * (anp.exp((1 - ρ) *vmk) * anp.exp((1 - ρ) * log_gk))
        foc = -(1-β)*anp.exp((1 - ρ)*log_cmk) * (1/(a-imk)) + β * ε * anp.exp((1 - ρ)* (log_gk + vmk)) * (1/(1+ ϕ_2 * imk))
        capital_growth = log_gk - ε * (ϕ_1 * anp.log(1. + ϕ_2 * imk) - α_k)
        output_constraint = log_cmk - anp.log(a - imk)

        return anp.array([recursive, foc, capital_growth, output_constraint])

    initial_guess = [-5, 0.0664766959220647, -3.86035793, 0.02330655]
    root = fsolve(f, initial_guess)

    vmk, imk, log_cmk, log_gk = root
    log_imk = anp.log(imk)

    #Assemble variables for equilibrium conditions
    rmv = log_gk
    log_y = anp.log(mu_2)

    if return_recursive:
        return anp.array([rmv, vmk, log_cmk, log_imk, log_gk, 0., log_y])
    else:
        return anp.array([log_cmk, log_imk, log_gk, 0., log_y])

2.4 Specifying Number of Variables #

We need to specify the number of jump variables, state and shock variables as an array (n_J, n_X, n_W).

2.5 Model Parameters #

Next, we need to specify all the model parameters using a Python tuple.

The first three parameters should be γ, β, ρ, respectively. The package is designed for recursive utlity.
Please use 1.00001 for γ and ρ as approximation for 1.0.
We suggest not using matrix when specifying the equilibrium conditions and parameters.

2.6 Additional Functions #

Finally, we specify how to compute the log-growth of consumption \(\log{C_{t+\epsilon}/C_t}\), log-growth of capital \(\log{G_{t+\epsilon}/G_t}\) and log-consumption-capital ratio \(\log{C_{t}/G_t}\) as a function of the defined variables Var_t and Var_tp1, which will be useful in computing the approximations for \(N_{t+\epsilon}\) and \(Q_t\). The log-growth of capital is already a state variable in this particular specification, but we list it here for more general purposes.

To compute log-growth of consumption, we simply use the decomposition:

\[ \log{C_{t+\epsilon}/C_t} = \log{C_{t+\epsilon}/K_{t+\epsilon}} + \log{K_{t+\epsilon}/K_t} - \log{C_t/K_t} \]

def gc_onecap_3d(Var_t, Var_tp1, W_tp1, q, *args):
    # Unpack model variables
    cmk_t, imk_t, gk_t, Z_t, Y_t = Var_t.ravel()
    cmk_tp1, imk_tp1, gk_tp1, Z_tp1, Y_tp1 = Var_tp1.ravel() 

    # Compute log consumption growth
    gc_tp1 = cmk_tp1 + gk_tp1 - cmk_t
    
    return gc_tp1

def gk_onecap_3d(Var_t, Var_tp1, W_tp1, q, *args):
    # Unpack model variables
    cmk_t, imk_t, gk_t, Z_t, Y_t = Var_t.ravel()
    cmk_tp1, imk_tp1, gk_tp1, Z_tp1, Y_tp1 = Var_tp1.ravel() 
    
    return gk_tp1

def cmk_onecap_3d(Var_t, Var_tp1, W_tp1, q, *args):
    # Unpack model variables
    cmk_t, imk_t, gk_t, Z_t, Y_t = Var_t.ravel()
    cmk_tp1, imk_tp1, gk_tp1, Z_tp1, Y_tp1 = Var_tp1.ravel() 
    
    return cmk_t

2.7 Computing ModelSol #

Now we are able to use steps 2.1 to 2.5 as inputs to the uncertain_expansion function. The solution is stored in a class named ModelSolution(Please refer to uncertainexpansion.ipynb for details.) under the Linear Quadratic Framework.

See src/uncertain_expansion.py for the source code of the expansion suite.

ModelSol['JX1_t'].coeffs

{'c': array([[-0.00000928],
        [ 0.00000359],
        [ 0.        ],
        [ 0.        ],
        [ 0.        ]]),
 'x': array([[ 0.        ,  0.17006158, -0.        ],
        [-0.        , -0.0658057 , -0.        ],
        [ 1.        ,  0.        ,  0.        ],
        [ 0.        ,  1.        ,  0.        ],
        [ 0.        ,  0.        ,  1.        ]])}

3 Outputs #

3.1 List of Outputs #

We now examine the contents of ModelSol, which contains the attributes listed below. Each approximation is stored in a class LinQuadVar, which contains the coefficients for \(X_{1,t}, X_{2,t}, X_{1,t}'X_{1,t}, W_{t+\epsilon}, W_{t+\epsilon}'W_{t+\epsilon}, X_{1,t}'W_{t+\epsilon}\) and the constant.

Input	Type	Description
`JXn_t`	LinQuadVar	Approximation of jump and state variables at time \(t\). Replace `n` with `0,1,2` to get the zeroth, first and second-order contribution. Omit `n` to get the full approximation. The variables are indexed in the order specified in Section 2.
`Jn_t`	LinQuadVar	Same as `JXn_t` but limited to jump variables.
`Xn_tp1`	LinQuadVar	Same as `JXn_t` but limited to state variables.
`JXn_t_tilde`	LinQuadVar	Same as `JXn_t` but using distorted measure. This variation is also available for `Jn_t` and `Xn_tp1`.
`util_sol`	dict	Dictionary containing solutions of the continuation values, including \(\mu^0, \Upsilon_0^2, \Upsilon_1^2,\) and \(\Upsilon_2^2\) etc.
`vmrn_tp1`	LinQuadVar	Approximation of continuation values \(\widehat{V}^1_{t+\epsilon}-\widehat{R}^1_t\) . Replace `n` with `0,1,2` to get the zeroth, first and second-order contribution. Omit `n` to get the full approximation.
`gcn_tp1`	LinQuadVar	Approximation of consumption growth \(\widehat{C}_{t+\epsilon}-\widehat{C}_t\) . Replace `n` with `0,1,2` to get the zeroth, first and second-order contribution. Omit `n` to get the full approximation.
`ss`	dict	Steady states for state and jump variables
`log_N_tilde`	LinQuadVar	Approximation for the log change of measure

For example, we can obtain the coefficients for the first-order contribution of \(\log{C_t/K_t}\) in the following way, noting that cmk_t was listed as the first jump variable when we specified the equilibrum conditions.

## Log consumption over capital approximation results, shown in the LinQuadVar coefficients form
ModelSol['JX1_t'][0].coeffs

{'c': array([[-0.00000928]]),
 'x': array([[ 0.        ,  0.17006158, -0.        ]])}

We can also display the full second-order approximation of \(\log{C_t/K_t}\) using the disp function, which renders a LinQuadVar object in LATEX form.

## Log consumption over capital approximation results, shown in the Latex analytical form
disp(ModelSol['JX_t'][0],'\\log\\frac{C_t}{K_t}') 

\[\begin{split}\displaystyle \log\frac{C_t}{K_t}=-3.661+\begin{bmatrix}2.473e-36&0.1701\end{bmatrix}X_t^1+\begin{bmatrix}2.576e-39&0.08503\end{bmatrix}X_t^2+X^{1T}_{t}\begin{bmatrix}-3.805e-51&-2.111e-38\\-7.732e-71&-1.608e-38\end{bmatrix}X^1_{t}\end{split}\]

Another example:

## Log capital growth process second order approximation results
disp(ModelSol['X2_tp1'][0],'\\log\\frac{K_{t+\epsilon}^2}{K_t^2}') 

\[\begin{split}\displaystyle \log\frac{K_{t+\epsilon}^2}{K_t^2}=-4.177e-05+\begin{bmatrix}8.061e-43&3.912e-08\end{bmatrix}X_t^1+\begin{bmatrix}-4.152e-38&0.03714\end{bmatrix}X_t^2+X^{1T}_{t}\begin{bmatrix}1.278e-52&-7.301e-39\\2.597e-72&-7.47e-39\end{bmatrix}X^1_{t}+X^{1T}_{t}\begin{bmatrix}0&0\\0&0\end{bmatrix}W_{t+1}\end{split}\]

3.3 Simulate Variables #

Given a series of shock processes, we can simulate the path of our state and jump variables using the ModelSolution.simulate method. Here, we simulate 400 periods of i.i.d standard multivariate normal shocks.

Ws = np.random.multivariate_normal(np.zeros(n_W),np.eye(n_W),size = 400)
JX_sim = ModelSol.simulate(Ws)

4 Using `LinQuadVar` in Computation #

In the previous section, we saw how to use uncertain_expansion to approximate variables and store their coefficients as LinQuadVar objects. In this section, we explore how to manipulate LinQuadVar objects for different uses.

To aid our examples, we first extract the steady states for the state evolution processes from the previous model solution:

See src/lin_quad.py for source code of LinQuadVar definition.

ModelSol['ss']

array([ -3.66035793,  -2.71090383,   0.01330655,   0.        , -11.97496092])

n_J, n_X, n_W = ModelSol['var_shape']
X0_tp1 = LinQuadVar({'c':np.array([[ModelSol['ss'][0]],[ModelSol['ss'][1]],[ModelSol['ss'][2]]])}, shape = (n_X, n_X, n_W))

4.1 `LinQuadVar` Operations #

We can sum multiple LinQuads together in two different ways. Here we demonstrate this with an example by summing the zeroth, first and second order contributions of our approximation for capital growth.

gk_tp1 = X0_tp1[0] + ModelSol['X1_tp1'][0]  + 0.5 * ModelSol['X2_tp1'][0] 
disp(gk_tp1,'\\log\\frac{K_{t+\epsilon}}{K_t}') 

\[\begin{split}\displaystyle \log\frac{K_{t+\epsilon}}{K_t}=-3.66+\begin{bmatrix}-4.152e-38&0.03714\end{bmatrix}X_t^1+\begin{bmatrix}0.007999&0.003478\end{bmatrix}W_{t+1}+\begin{bmatrix}-2.076e-38&0.01857\end{bmatrix}X_t^2+X^{1T}_{t}\begin{bmatrix}6.389e-53&-3.651e-39\\1.298e-72&-3.735e-39\end{bmatrix}X^1_{t}+X^{1T}_{t}\begin{bmatrix}0&0\\0&0\end{bmatrix}W_{t+1}\end{split}\]

In the next example, we sum together the contributions for both capital growth and technology:

lq_sum([X0_tp1, ModelSol['X1_tp1'], 0.5 * ModelSol['X2_tp1']]).coeffs

{'x2': array([[-0.       ,  0.0185721,  0.       ],
        [ 0.       ,  0.472    ,  0.       ],
        [ 0.       ,  0.       ,  0.403    ]]),
 'x': array([[-0.        ,  0.03714422,  0.00000001],
        [ 0.        ,  0.944     ,  0.        ],
        [ 0.        ,  0.        ,  0.806     ]]),
 'c': array([[-3.66037866],
        [-2.71090383],
        [-0.07821726]]),
 'xw': array([[ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,  0.        ,  0.00399962,  0.00173897,  0.        ],
        [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,  0.        ,  0.        ,  0.02478025,  0.        ],
        [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,  0.        ,  0.        ,  0.        , -0.21392032]]),
 'w': array([[0.00799924, 0.00347793, 0.        ],
        [0.        , 0.04956051, 0.        ],
        [0.        , 0.        , 0.42784065]]),
 'xx': array([[ 0.        , -0.        ,  0.        , -0.        , -0.00021926,  0.        ,  0.        ,  0.        ,  0.        ],
        [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,  0.        ,  0.        ,  0.        ,  0.        ],
        [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,  0.        ,  0.        ,  0.        ,  0.097     ]])}

4.2 `LinQuadVar` Split and Concat #

split breaks multiple dimensional LinQuad into one-dimensional LinQuads, while concat does the inverse.

gk_tp1, Z_tp1, Y_tp1 = ModelSol['X1_tp1'].split()
concat([gk_tp1, Z_tp1, Y_tp1])

<lin_quad.LinQuadVar at 0x314574c50>

4.3 Evaluate a `LinQuadVar` #

We can evaluate a LinQuad at specific state \((X_{t},W_{t+\epsilon})\) in time. As an example, we evaluate all 5 variables under steady state with a multivariate random normal shock.

x1 = np.zeros([n_X ,1])
x2 = np.zeros([n_X ,1])
w = np.random.multivariate_normal(np.zeros(n_W),np.eye(n_W),size = 1).T
ModelSol['JX_tp1'](*(x1,x2,w))

array([[ -3.66764679],
       [ -2.70809391],
       [  0.00945342],
       [ -0.03677232],
       [-12.41358991]])

4.4 Next period expression for `LinQuadVar` #

ModelSol allows us to express a jump variable \(J_t\) as a function of \(t\) state and shock variables. Suppose we would like to compute its next period expression \(J_{t+\epsilon}\) with shocks. The function next_period expresses \(J_{t+\epsilon}\) in terms of \(t\) state variables and \(t+\epsilon\) shock variables. For example, we can express the \(t+\epsilon\) expression for the first-order contribution to consumption over capital as:

cmk1_tp1 = next_period(ModelSol['J1_t'][0], ModelSol['X1_tp1'])
disp(cmk1_tp1, '\\log\\frac{C_{t+\epsilon}^1}{K_{t+\epsilon}^1}') 

\[\displaystyle \log\frac{C_{t+\epsilon}^1}{K_{t+\epsilon}^1}=-9.281e-06+\begin{bmatrix}-1.027e-73&0.1605\end{bmatrix}X_t^1+\begin{bmatrix}1.978e-38&0.008428\end{bmatrix}W_{t+1}\]

cmk2_tp1 = next_period(ModelSol['J2_t'][0], ModelSol['X1_tp1'], ModelSol['X2_tp1'])
disp(cmk2_tp1, '\\log\\frac{C_{t+\epsilon}^2}{K_{t+\epsilon}^2}') 

\[\begin{split}\displaystyle \log\frac{C_{t+\epsilon}^2}{K_{t+\epsilon}^2}=-0.002043+\begin{bmatrix}-1.644e-79&1.181e-06\end{bmatrix}X_t^1+\begin{bmatrix}3.247e-44&6.201e-08\end{bmatrix}W_{t+1}+\begin{bmatrix}-2.14e-76&0.1605\end{bmatrix}X_t^2+X^{1T}_{t}\begin{bmatrix}6.584e-91&1.618e-75\\5.189e-108&1.222e-75\end{bmatrix}X^1_{t}+X^{1T}_{t}\begin{bmatrix}5.055e-90&1.531e-76\\2.747e-108&-5.617e-40\end{bmatrix}W_{t+1}+W_{t+1}^{T}\begin{bmatrix}-4.869e-55&-1.674e-41\\-5.292e-73&-1.275e-41\end{bmatrix}W_{t+1}\end{split}\]

4.6 Compute the Expectation of time \(t+\epsilon\) `LinQuadVar` #

Suppose the distribution of shocks has a constant mean and variance (not state dependent). Then, we can use the E function to compute the expectation of a time \(t+\epsilon\) LinQuadVar as follows:

E_w = ModelSol['util_sol']['μ_0']
cov_w = np.eye(n_W)
E_ww = cal_E_ww(E_w, cov_w)
E_cmk2_tp1 = E(cmk2_tp1, E_w, E_ww)
disp(E_cmk2_tp1, '\mathbb{E}[\\log\\frac{C_{t+\epsilon}^2}{K_{t+\epsilon}^2}|\mathfrak{F_t}]')

\[\begin{split}\displaystyle \mathbb{E}[\log\frac{C_{t+\epsilon}^2}{K_{t+\epsilon}^2}|\mathfrak{F_t}]=-0.00206+\begin{bmatrix}-1.691e-79&1.2e-06\end{bmatrix}X_t^1+\begin{bmatrix}-2.14e-76&0.1605\end{bmatrix}X_t^2+X^{1T}_{t}\begin{bmatrix}6.584e-91&1.618e-75\\5.189e-108&1.222e-75\end{bmatrix}X^1_{t}\end{split}\]

Suppose the distribution of shock has a state-dependent mean and variance (implied by \(\tilde{N}_{t+\epsilon}\) shown in the notes), we can use E_N_tp1 and N_tilde_measure to compute the expectation of time \(t+\epsilon\) LinQuadVar.

N_cm = N_tilde_measure(ModelSol['util_sol']['log_N_tilde'],(1,n_X,n_W))
E_cmk2_tp1_tilde = E_N_tp1(cmk2_tp1, N_cm)
disp(E_cmk2_tp1_tilde, '\mathbb{\\tilde{E}}[\\log\\frac{C_{t+\epsilon}^2}{K_{t+\epsilon}^2}|\mathfrak{F_t}]')

\[\begin{split}\displaystyle \mathbb{\tilde{E}}[\log\frac{C_{t+\epsilon}^2}{K_{t+\epsilon}^2}|\mathfrak{F_t}]=-0.00206+\begin{bmatrix}-1.691e-79&1.2e-06\end{bmatrix}X_t^1+\begin{bmatrix}-2.14e-76&0.1605\end{bmatrix}X_t^2+X^{1T}_{t}\begin{bmatrix}6.584e-91&1.618e-75\\-4.268e-80&1.222e-75\end{bmatrix}X^1_{t}\end{split}\]

Uncertain Expansion Example Notebook

Contents

Uncertain Expansion Example Notebook#

1 Model Setup#