{
"cells": [
{
"cell_type": "markdown",
"id": "f5e7b1d0",
"metadata": {
"id": "f5e7b1d0"
},
"source": [
"
"
]
},
{
"cell_type": "markdown",
"id": "7735eea6",
"metadata": {},
"source": [
"# Uncertain Expansion Example Notebook"
]
},
{
"cell_type": "markdown",
"id": "ef624dcb",
"metadata": {
"id": "ef624dcb"
},
"source": [
"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).\n",
"\n",
"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.\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "6f590ad3-c8d4-4fba-93f4-bf14f4f3364b",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"```{contents}\n",
":depth: 1\n",
"```\n",
"\n",
"$\\newcommand{\\eqdef}{\\stackrel{\\text{def}}{=}}$"
]
},
{
"cell_type": "markdown",
"id": "6f6c2319",
"metadata": {},
"source": [
"# 1 Model Setup"
]
},
{
"cell_type": "markdown",
"id": "1141ee3e",
"metadata": {},
"source": [
"We employ a discrete-time approximation of the single-capital model in Section 4.4 of {cite}`HansenKhorramiTourre:2024`, as outlined in Section 11.7.1.\n",
"\n",
"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}$.\n",
"\n",
"The law of motions for the three states are:\n",
"```{math}\n",
"\\begin{align}\n",
"{{\\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) \n",
" + \\sqrt{\\epsilon} \\exp\\left(\\frac 1 2 { Z}_{2,t} \\right) {\\sigma}_k W_{t+\\epsilon} \n",
" \\end{align} \n",
"```\n",
"\n",
"$$\n",
"\\begin{align}\n",
"{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) \n",
"+ \\sqrt{\\epsilon} \\exp\\left( - {\\frac 1 2} {Z}_{2,t} \\right) \\sigma_2 W_{t+\\epsilon}. \\end{align} \n",
"$$\n",
"\n",
"$$\n",
"\\begin{align}\n",
"{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) \n",
"+ \\sqrt{\\epsilon} \\exp\\left( - {\\frac 1 2} {Z}_{2,t} \\right) \\sigma_2 W_{t+\\epsilon}.\n",
"\\end{align} \n",
"$$\n",
"\n",
"We can solve for the jump variables using the first-order condition for $D_t$ and the output constraint:\n",
"\n",
"$$\n",
"{\\widehat C}_t = \\log \\left( \\alpha - D_t \\right) + {\\widehat K}_t \n",
"$$\n",
"$$\n",
"\\begin{align*}\n",
"0 &= - (1-\\beta) \\exp\\left[(1- \\rho) \\left({\\widehat C}_t - {\\widehat G}_t \\right)\\right] \\left( \\frac{1}{\\alpha - D_t} \\right) \n",
"\\\\&+ \\beta \\exp\\left[(1 - \\rho) \\left({\\widehat R}_t - {\\widehat G}_t \\right)\\right] \\left(\\frac{1}{1 + \\zeta D_t} \\right) \n",
"\\end{align*}\n",
"$$\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 1.2 Rewriting First-Order Conditions \n",
"\n",
"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\n",
"\n",
"$$\n",
"\\begin{align*}\n",
"Q_t \\eqdef & \\hspace{.2cm} \\beta \\exp\\left[(1-\\rho) \\left( {\\widehat R}_t - {\\widehat V}_t \\right) \\right] \\cr \\cr\n",
"P_t \\eqdef & \\hspace{.2cm} (1-\\beta) \\exp\\left[(\\rho-1) \\left( {\\widehat V}_t - {\\widehat G}_t \\right) \\right] \\cr \\cr\n",
"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 \n",
"\\psi_{g'}^x (D_t, X_t, W_{t+\\epsilon})' & 1 \\end{bmatrix} \\begin{bmatrix} MX_{t+\\epsilon} \\cr MG_{t+\\epsilon} \\end{bmatrix} \n",
"\\ \\cr \\cr\n",
"L_t \\eqdef & \\hspace{.2cm} \\begin{bmatrix} \\exp \\left[ (1-\\rho) \\left({\\widehat C}_t - {\\widehat G}_t \\right) \\right] \n",
" \\begin{bmatrix} \\kappa_d(D_t,X_t) \\cr \\kappa_x(D_t, X_t) \\cr 1 \\end{bmatrix} \\end{bmatrix} \\cr\n",
"M_t \\eqdef & \\hspace{.2cm} \\begin{bmatrix} 0 \\cr MG_t \\cr MX_t \\end{bmatrix} .\n",
"\\end{align*}\n",
"$$\n",
"\n",
"\n",
"\n",
"By noting that $MG_t=1$, we can express the model above simply as:\n",
"\n",
"$$\n",
"H_{t+\\epsilon} = \\begin{bmatrix} 0 \\\\\n",
"\\frac \\epsilon {1+\\zeta D_t} \n",
"\\end{bmatrix}\n",
"$$\n",
"\n",
"$$\n",
"L_{t} = \\begin{bmatrix}\n",
"0 \\\\ \\exp\\left[(1-\\rho)(\\widehat{C}_t-\\widehat{G}_t)\\right] \\frac 1 {\\alpha-D_t}\n",
"\\end{bmatrix}\n",
"$$\n",
"\n",
"$$\n",
"U_t = \\begin{bmatrix} \\exp{(\\widehat{C}_t - \\widehat{K}_t)} + D_t - \\alpha \\\\\n",
"0\n",
"\\end{bmatrix}\n",
"$$\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"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. \n",
"\n",
"\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"id": "f5e59370",
"metadata": {},
"source": [
"
\n",
"
\n",
"
\n",
"\n",
"# 2 Inputs \n",
"\n",
"The `Expansion Suite` uses the function `uncertain_expansion` to approximate a solution to the above system locally. It takes the following inputs:\n",
"\n",
"```{list-table}\n",
":header-rows: 1\n",
"\n",
"* - Input\n",
" - Type\n",
" - Description\n",
"* - `eq`\n",
" - *callable*\n",
" - Function which takes state and jump variables and outputs error of equilibrium conditions, i.e. RHS of equation (1)\n",
"* - `ss`\n",
" - *callable*\n",
" - Function which takes model parameters and outputs steady-state values of state and jump variables\n",
"* - `var_shape`\n",
" - *tuple of int*\n",
" - Triple of integers $(n_J,n_X,n_W)$ where $n_J$ is the number of jump variables, $n_X$ is the number of state variables, and $n_W$ is the number of shocks.\n",
"* - `args`\n",
" - *tuple of float or ndarray*\n",
" - Array of preference and model parameters\n",
"* - `gc`\n",
" - *callable*\n",
" - Function which takes in variables and outputs log-growth of consumption\n",
"* - `args`\n",
" - tuple of float/ndarray\n",
" - Preference and model parameters\n",
"* - `gc`\n",
" - *callable*\n",
" - Function which takes in variables and outputs log-growth of consumption\n",
"\n",
"* - `init_util`\n",
" - *dict, or None*\n",
" - Dictionary for initialization of the expansion coefficients $\\mu^0, \\Upsilon_0^2, \\Upsilon_1^2,$ and $\\Upsilon_2^2$. The keys are: `mu_0`, `Upsilon_0`, `Upsilon_1`, `Upsilon_2`. If *None*, zero matrices are used\n",
"* - `iter_tol`\n",
" - *float*\n",
" - Iteration stops when the maximum absolute difference between the approximated variables in the current and previous iteration is below `iter_tol`\n",
"* - `max_iter`\n",
" - *int*\n",
" - Maximum number of iterations\n",
"```\n",
"\n",
"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. \n",
"\n",
"We will now walk through the computation using the example above.\n",
"
\n",
"\n",
"## 2.1 Libraries \n",
"\n",
"Begin by installing the following libraries and downloading `RiskUncertaintyValue`, which contains the functions required to solve the model:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "f0d30383",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"fatal: destination path 'RiskUncertaintyValue' already exists and is not an empty directory.\n"
]
},
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import os\n",
"import sys\n",
"workdir = os.getcwd()\n",
"!git clone https://github.com/lphansen/RiskUncertaintyValue \n",
"workdir = os.getcwd() + '/RiskUncertaintyValue' \n",
"sys.path.insert(0, workdir+'/src') \n",
"import numpy as np\n",
"import autograd.numpy as anp\n",
"from scipy import optimize\n",
"np.set_printoptions(suppress=True)\n",
"np.set_printoptions(linewidth=200)\n",
"from IPython.display import display, HTML\n",
"from BY_example_sol import disp\n",
"display(HTML(\"\"))\n",
"import warnings\n",
"warnings.filterwarnings(\"ignore\")\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from lin_quad import LinQuadVar\n",
"from lin_quad_util import E, concat, next_period, cal_E_ww, lq_sum, N_tilde_measure, E_N_tp1\n",
"from uncertain_expansion import uncertain_expansion\n",
"np.set_printoptions(suppress=True)\n",
"import pickle\n",
"from scipy.optimize import fsolve"
]
},
{
"cell_type": "markdown",
"id": "06a5ca90-4eec-48d3-862a-4a00d766b58e",
"metadata": {},
"source": [
"## 2.2 Equilibrium Condition Function \n",
"\n",
"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:\n",
"\n",
"```{list-table}\n",
":header-rows: 1\n",
"\n",
"* - Variable\n",
" - Type\n",
" - Description\n",
" - Corresponds to\n",
"* - `Var_t`\n",
" - array-like\n",
" - Vector of jump and state variables in the current period\n",
" - $(J_t, X_t)$\n",
"* - `Var_tp1`\n",
" - array-like\n",
" - Vector of jump and state variables in the next period\n",
" - $(J_{t+\\epsilon}, X_{t+\\epsilon})$\n",
"* - `W_tp1`\n",
" - array-like\n",
" - Vector of shock variables in the next period\n",
" - $(W_{t+\\epsilon})$\n",
"* - `q`\n",
" - float\n",
" - Perturbation parameter\n",
" - $q$\n",
"* - `mode`\n",
" - string\n",
" - 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\n",
" - \n",
"* - `recursive_ss`\n",
" - array-like\n",
" - Steady-state values of $N$ and $Q$\n",
" - \n",
"* - `args`\n",
" - tuple of float/ndarray\n",
" - Preference and model parameters\n",
" - \n",
"```\n"
]
},
{
"cell_type": "markdown",
"id": "ba75562f",
"metadata": {
"id": "ba75562f"
},
"source": [
"\n",
"While writing out `eq`, ensure that: \n",
"\n",
"\n",
"- State variables follow jump variables.\n",
"- The number of jump variables equals the number of forward-looking equilibrium conditions. \n",
"- The number of state variables equals to the number of state evolution equations.\n",
"- All values are expressed as type `float`. Otherwise, this may cause type errors.\n",
"\n",
"\n",
"For the example model, we write:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "6ae5d2d3",
"metadata": {
"id": "6ae5d2d3"
},
"outputs": [],
"source": [
"def eq_onecap_3d(Var_t, Var_tp1, W_tp1, q, mode, recursive_ss, *args):\n",
" \n",
" # Unpack model parameters\n",
" γ, β, ρ, 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\n",
" w1_tp1, w2_tp1, w3_tp1 = W_tp1.ravel()\n",
"\n",
" # Unpack model variables\n",
" log_cmk_t, log_imk_t, gk_t, Z_t, Y_t = Var_t.ravel()\n",
" log_cmk_tp1, log_imk_tp1, gk_tp1, Z_tp1, Y_tp1 = Var_tp1.ravel() \n",
"\n",
" imk_t = anp.exp(log_imk_t)\n",
" cmk_t = anp.exp(log_cmk_t)\n",
"\n",
" q_t = β*anp.exp((1.-ρ)*recursive_ss[0])\n",
" p_t = (1-β)*anp.exp((ρ-1.)*recursive_ss[1])\n",
"\n",
" # State evolution processes\n",
" 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)\n",
" # Equation for states at time t\n",
" phi_1 = gk_tp1 - dgk\n",
" 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)\n",
" phi_3 = Y_tp1 - Y_t + ε * beta_2 * (1 - mu_2 * anp.exp(-Y_t)) \\\n",
" + q**2 * 0.5 * (sigma_y1**2 + sigma_y2**2 + sigma_y3**2) * anp.exp(-Y_t) * ε \\\n",
" - anp.exp(-0.5 * Y_t) * (sigma_y1*w1_tp1+sigma_y2*w2_tp1+sigma_y3*w3_tp1) * anp.sqrt(ε)\n",
"\n",
" ## Output constraint\n",
" qnh_1 = 0.\n",
" pl_1 = 0.\n",
" u_1 = -a + cmk_t + imk_t\n",
"\n",
" ## Investment FOC\n",
" qnh_2 = ε/(1.+ϕ_2*imk_t) \n",
" pl_2 = - anp.exp((1.-ρ)*log_cmk_t) * 1/(a-imk_t)\n",
" u_2 = 0.\n",
"\n",
" if mode == 'H':\n",
" return anp.array([qnh_1, qnh_2]) \n",
" \n",
" if mode == 'L':\n",
" return anp.array([pl_1, pl_2]) \n",
"\n",
" return anp.array([\n",
" q_t*qnh_1 - u_1,\n",
" q_t*qnh_2 + p_t* pl_2 - u_2,\n",
" phi_1, phi_2, phi_3])\n"
]
},
{
"cell_type": "markdown",
"id": "d554abc5",
"metadata": {
"id": "d554abc5"
},
"source": [
"## 2.3 Steady State Function \n",
"\n",
"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. \n",
"\n",
""
]
},
{
"cell_type": "code",
"execution_count": 44,
"id": "e0d871c4",
"metadata": {
"id": "e0d871c4"
},
"outputs": [],
"source": [
"def ss_onecap_3d(*args, return_recursive=False):\n",
" γ, β, ρ, 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\n",
"\n",
" def f(x):\n",
" vmk, imk, log_cmk, log_gk = x\n",
"\n",
" recursive = anp.exp((1 - ρ) * vmk) - (1 - β) * anp.exp((1 - ρ) * log_cmk) - β * (anp.exp((1 - ρ) *vmk) * anp.exp((1 - ρ) * log_gk))\n",
" foc = -(1-β)*anp.exp((1 - ρ)*log_cmk) * (1/(a-imk)) + β * ε * anp.exp((1 - ρ)* (log_gk + vmk)) * (1/(1+ ϕ_2 * imk))\n",
" capital_growth = log_gk - ε * (ϕ_1 * anp.log(1. + ϕ_2 * imk) - α_k)\n",
" output_constraint = log_cmk - anp.log(a - imk)\n",
"\n",
" return anp.array([recursive, foc, capital_growth, output_constraint])\n",
"\n",
" initial_guess = [-5, 0.0664766959220647, -3.86035793, 0.02330655]\n",
" root = fsolve(f, initial_guess)\n",
"\n",
" vmk, imk, log_cmk, log_gk = root\n",
" log_imk = anp.log(imk)\n",
"\n",
" #Assemble variables for equilibrium conditions\n",
" rmv = log_gk\n",
" log_y = anp.log(mu_2)\n",
"\n",
" if return_recursive:\n",
" return anp.array([rmv, vmk, log_cmk, log_imk, log_gk, 0., log_y])\n",
" else:\n",
" return anp.array([log_cmk, log_imk, log_gk, 0., log_y])\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "aadc21b2",
"metadata": {
"id": "aadc21b2"
},
"source": [
"## 2.4 Specifying Number of Variables \n",
"\n",
"We need to specify the number of jump variables, state and shock variables as an array `(n_J, n_X, n_W)`."
]
},
{
"cell_type": "code",
"execution_count": 45,
"id": "6f4f50b9",
"metadata": {
"id": "6f4f50b9",
"tags": [
"hide-cell"
]
},
"outputs": [],
"source": [
"var_shape = (2, 3, 3)"
]
},
{
"cell_type": "markdown",
"id": "715b17d5",
"metadata": {
"id": "715b17d5"
},
"source": [
"## 2.5 Model Parameters \n",
"\n",
"Next, we need to specify all the model parameters using a Python tuple. \n",
"\n",
"- The first three parameters should be γ, β, ρ, respectively. The package is designed for recursive utlity.\n",
"- Please use 1.00001 for γ and ρ as approximation for 1.0.\n",
"- We suggest not using matrix when specifying the equilibrium conditions and parameters.\n",
"\n",
""
]
},
{
"cell_type": "code",
"execution_count": 92,
"id": "43a4ac46",
"metadata": {
"id": "43a4ac46",
"tags": [
"hide-cell"
]
},
"outputs": [],
"source": [
"#Params\n",
"sigma_k1 = 0.92 * anp.sqrt(12)\n",
"sigma_k2 = 0.4 * anp.sqrt(12)\n",
"sigma_k3 = 0.0\n",
"sigma_z1 = 0.0\n",
"sigma_z2 = 5.7 * anp.sqrt(12)\n",
"sigma_z3 = 0.0\n",
"sigma_y1 = 0.0\n",
"sigma_y2 = 0.0\n",
"sigma_y3 = 0.00031 * anp.sqrt(12)\n",
"\n",
"\n",
"delta = 0.01\n",
"alpha = 0.0922\n",
"\n",
"a = alpha\n",
"ϕ_1 = 1/8\n",
"ϕ_2 = 8.\n",
"beta_k = 0.04\n",
"beta_z = 0.056\n",
"beta_2 = 0.194\n",
"ε = 1.0\n",
"mu_2 = 6.3 * (10**(-6))\n",
"\n",
"α_k = 0.04\n",
"γ = 8.0\n",
"ρ = 1.001\n",
"β = anp.exp(-ε * delta)\n",
"\n",
"args = γ, β, ρ, 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, ε \n"
]
},
{
"cell_type": "markdown",
"id": "238e1c2a",
"metadata": {
"id": "238e1c2a"
},
"source": [
"## 2.6 Additional Functions \n",
"\n",
"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. \n",
"\n",
"To compute log-growth of consumption, we simply use the decomposition:\n",
"\n",
"$$\n",
"\\log{C_{t+\\epsilon}/C_t} = \\log{C_{t+\\epsilon}/K_{t+\\epsilon}} + \\log{K_{t+\\epsilon}/K_t} - \\log{C_t/K_t}\n",
"$$\n",
"\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 93,
"id": "961ab479",
"metadata": {},
"outputs": [],
"source": [
"def gc_onecap_3d(Var_t, Var_tp1, W_tp1, q, *args):\n",
" # Unpack model variables\n",
" cmk_t, imk_t, gk_t, Z_t, Y_t = Var_t.ravel()\n",
" cmk_tp1, imk_tp1, gk_tp1, Z_tp1, Y_tp1 = Var_tp1.ravel() \n",
"\n",
" # Compute log consumption growth\n",
" gc_tp1 = cmk_tp1 + gk_tp1 - cmk_t\n",
" \n",
" return gc_tp1\n",
"\n",
"def gk_onecap_3d(Var_t, Var_tp1, W_tp1, q, *args):\n",
" # Unpack model variables\n",
" cmk_t, imk_t, gk_t, Z_t, Y_t = Var_t.ravel()\n",
" cmk_tp1, imk_tp1, gk_tp1, Z_tp1, Y_tp1 = Var_tp1.ravel() \n",
" \n",
" return gk_tp1\n",
"\n",
"def cmk_onecap_3d(Var_t, Var_tp1, W_tp1, q, *args):\n",
" # Unpack model variables\n",
" cmk_t, imk_t, gk_t, Z_t, Y_t = Var_t.ravel()\n",
" cmk_tp1, imk_tp1, gk_tp1, Z_tp1, Y_tp1 = Var_tp1.ravel() \n",
" \n",
" return cmk_t\n"
]
},
{
"cell_type": "markdown",
"id": "020b982b",
"metadata": {
"id": "020b982b"
},
"source": [
"## 2.7 Computing ModelSol \n",
"\n",
"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. \n",
"\n",
"See [src/uncertain_expansion.py](https://github.com/lphansen/RiskUncertaintyValue/blob/main/src/uncertain_expansion.py) for the source code of the expansion suite.\n"
]
},
{
"cell_type": "code",
"execution_count": 96,
"id": "a758913c",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"{'c': array([[-0.00000928],\n",
" [ 0.00000359],\n",
" [ 0. ],\n",
" [ 0. ],\n",
" [ 0. ]]),\n",
" 'x': array([[ 0. , 0.17006158, -0. ],\n",
" [-0. , -0.0658057 , -0. ],\n",
" [ 1. , 0. , 0. ],\n",
" [ 0. , 1. , 0. ],\n",
" [ 0. , 0. , 1. ]])}"
]
},
"execution_count": 96,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ModelSol['JX1_t'].coeffs"
]
},
{
"cell_type": "code",
"execution_count": 48,
"id": "217c9448",
"metadata": {
"id": "217c9448",
"outputId": "aab57b59-6dce-4452-dfde-a6ee003cdbd7",
"tags": [
"hide-cell"
]
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Iteration 1: error = 0.170925932\n",
"Iteration 2: error = 0.00105362328\n",
"Iteration 3: error = 4.51132298e-06\n",
"Iteration 4: error = 2.31767169e-08\n",
"Iteration 5: error = 1.19069643e-10\n"
]
}
],
"source": [
"# Solve the One-Capital Adjustment Model\n",
"eq = eq_onecap_3d\n",
"ss = ss_onecap_3d\n",
"gc = gc_onecap_3d\n",
"gk = gk_onecap_3d\n",
"cmk_app = cmk_onecap_3d\n",
"approach = '1' # See Exploring Recursive Utility Appendix for difference about approach '1' and '2'\n",
"ModelSol = uncertain_expansion(eq, ss, var_shape, args, gc, gk, cmk_app, approach = '1')"
]
},
{
"cell_type": "markdown",
"id": "961bdc55",
"metadata": {
"id": "961bdc55"
},
"source": [
"
\n",
"
\n",
"
\n",
"\n",
"# 3 Outputs\n",
"\n",
"## 3.1 List of Outputs \n",
"\n",
"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.\n",
"\n",
"\n",
"```{list-table}\n",
":header-rows: 1\n",
"\n",
"* - Input\n",
" - Type\n",
" - Description\n",
"* - `JXn_t`\n",
" - *LinQuadVar*\n",
" - 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. \n",
"* - `Jn_t`\n",
" - *LinQuadVar*\n",
" - Same as `JXn_t` but limited to jump variables.\n",
"* - `Xn_tp1`\n",
" - *LinQuadVar*\n",
" - Same as `JXn_t` but limited to state variables.\n",
"* - `JXn_t_tilde`\n",
" - *LinQuadVar*\n",
" - Same as `JXn_t` but using distorted measure. This variation is also available for `Jn_t` and `Xn_tp1`.\n",
"* - `util_sol`\n",
" - *dict*\n",
" - Dictionary containing solutions of the continuation values, including $\\mu^0, \\Upsilon_0^2, \\Upsilon_1^2,$ and $\\Upsilon_2^2$ etc.\n",
"* - `vmrn_tp1`\n",
" - *LinQuadVar*\n",
" - 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. \n",
"* - `gcn_tp1`\n",
" - *LinQuadVar*\n",
" - 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. \n",
"* - `ss`\n",
" - *dict*\n",
" - Steady states for state and jump variables\n",
"* - `log_N_tilde`\n",
" - *LinQuadVar*\n",
" - Approximation for the log change of measure\n",
"```\n",
"\n",
"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.\n"
]
},
{
"cell_type": "code",
"execution_count": 51,
"id": "3fa57ffa",
"metadata": {
"id": "3fa57ffa",
"outputId": "1cdb816b-4ed7-4c19-e009-60fdcf049004",
"tags": []
},
"outputs": [
{
"data": {
"text/plain": [
"{'c': array([[-0.00000928]]),\n",
" 'x': array([[ 0. , 0.17006158, -0. ]])}"
]
},
"execution_count": 51,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"## Log consumption over capital approximation results, shown in the LinQuadVar coefficients form\n",
"ModelSol['JX1_t'][0].coeffs"
]
},
{
"cell_type": "markdown",
"id": "6b4f132c",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": 52,
"id": "a21c4b7d",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\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}$"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"## Log consumption over capital approximation results, shown in the Latex analytical form\n",
"disp(ModelSol['JX_t'][0],'\\\\log\\\\frac{C_t}{K_t}') "
]
},
{
"cell_type": "markdown",
"id": "2c2d4266",
"metadata": {},
"source": [
"Another example:"
]
},
{
"cell_type": "code",
"execution_count": 53,
"id": "d5b18ec2",
"metadata": {
"id": "d5b18ec2",
"outputId": "1b2ebb60-d294-4281-e573-a76dcbc55b4d"
},
"outputs": [
{
"data": {
"text/latex": [
"$\\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}$"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"## Log capital growth process second order approximation results\n",
"disp(ModelSol['X2_tp1'][0],'\\\\log\\\\frac{K_{t+\\epsilon}^2}{K_t^2}') "
]
},
{
"cell_type": "markdown",
"id": "0190ee8e",
"metadata": {
"id": "0190ee8e"
},
"source": [
"\n",
"## 3.2 Approximate Related Variables \n",
"\n",
"We can use the method `ModelSol.approximate` to approximate variables which can be expressed as functions of the state and jump variables. For example, suppose we want to approximate log investment growth $\\log{I_{t+\\epsilon}/I_t}$. We can express this as:\n",
"\n",
"$$\n",
"\\log{I_{t+\\epsilon}/I_t} = \\log{I_{t+\\epsilon}/K_{t+\\epsilon}} + \\log{K_{t+\\epsilon}/K_t} - \\log{I_t/K_t}\n",
"$$\n",
"\n",
"We use this expression to write the function `gi_tp1_approx`, which is used as an input in `ModelSol.approximate` below."
]
},
{
"cell_type": "code",
"execution_count": 54,
"id": "145a47f6",
"metadata": {
"id": "145a47f6",
"outputId": "ae3e6f0d-39c8-47a6-b2b3-934eaa6cffaa",
"tags": [
"hide-cell"
]
},
"outputs": [
{
"data": {
"text/plain": [
"{'w': array([[-0. , -0.00326136, 0.42784065]]),\n",
" 'x': array([[0. , 0.00368512, 0.806 ]])}"
]
},
"execution_count": 54,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"def gi_tp1_approx(Var_t, Var_tp1, W_tp1, q, *args):\n",
"\n",
" # Parameters for the model\n",
" γ, β, ρ, 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\n",
"\n",
" # Variables: (q_t, q_tp1 is excluded when using the method in `ModelSolution`)\n",
" cmk_t, imk_t, mkmc_t, Z_t, gk_t = Var_t.ravel()\n",
" cmk_tp1, imk_tp1, mkmc_tp1, Z_tp1, gk_tp1 = Var_tp1.ravel()\n",
" \n",
" gi_tp1 = imk_tp1 + gk_tp1 -imk_t\n",
" \n",
" return gi_tp1\n",
"\n",
"n_J, n_X, n_W = var_shape\n",
"gi_tp1_list = ModelSol.approximate(gi_tp1_approx, args)\n",
"gi_tp1 = {'gi_tp1':gi_tp1_list[0],\\\n",
" 'gi0_tp1':gi_tp1_list[1],\\\n",
" 'gi1_tp1':gi_tp1_list[2],\\\n",
" 'gi2_tp1':gi_tp1_list[3]}\n",
"gi_tp1['gi1_tp1'].coeffs"
]
},
{
"cell_type": "markdown",
"id": "bbd28ede",
"metadata": {
"id": "bbd28ede"
},
"source": [
"
\n",
"\n",
"## 3.3 Simulate Variables \n",
"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."
]
},
{
"cell_type": "code",
"execution_count": 55,
"id": "58417c59",
"metadata": {
"id": "58417c59"
},
"outputs": [],
"source": [
"Ws = np.random.multivariate_normal(np.zeros(n_W),np.eye(n_W),size = 400)\n",
"JX_sim = ModelSol.simulate(Ws)"
]
},
{
"cell_type": "markdown",
"id": "eab36438",
"metadata": {
"id": "eab36438"
},
"source": [
"
\n",
"
\n",
"\n",
"# 4 Using `LinQuadVar` in Computation \n",
"\n",
"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.\n",
"\n",
"To aid our examples, we first extract the steady states for the state evolution processes from the previous model solution:\n",
"\n",
"See [src/lin_quad.py](https://github.com/lphansen/RiskUncertaintyValue/blob/main/src/lin_quad.py) for source code of `LinQuadVar` definition."
]
},
{
"cell_type": "code",
"execution_count": 56,
"id": "9698f3e2",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([ -3.66035793, -2.71090383, 0.01330655, 0. , -11.97496092])"
]
},
"execution_count": 56,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ModelSol['ss']"
]
},
{
"cell_type": "code",
"execution_count": 57,
"id": "b5d7ee6b",
"metadata": {
"id": "b5d7ee6b",
"outputId": "882669aa-3be7-4c20-dd54-1273a7c1826d"
},
"outputs": [],
"source": [
"n_J, n_X, n_W = ModelSol['var_shape']\n",
"X0_tp1 = LinQuadVar({'c':np.array([[ModelSol['ss'][0]],[ModelSol['ss'][1]],[ModelSol['ss'][2]]])}, shape = (n_X, n_X, n_W))"
]
},
{
"cell_type": "markdown",
"id": "a5270666",
"metadata": {
"id": "a5270666"
},
"source": [
"## 4.1 `LinQuadVar` Operations \n",
"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. "
]
},
{
"cell_type": "code",
"execution_count": 58,
"id": "205dd148",
"metadata": {
"id": "205dd148",
"outputId": "75bfa88f-3b99-4535-cf0d-121591eae233"
},
"outputs": [
{
"data": {
"text/latex": [
"$\\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}$"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"gk_tp1 = X0_tp1[0] + ModelSol['X1_tp1'][0] + 0.5 * ModelSol['X2_tp1'][0] \n",
"disp(gk_tp1,'\\\\log\\\\frac{K_{t+\\epsilon}}{K_t}') "
]
},
{
"cell_type": "markdown",
"id": "f901b7a3",
"metadata": {},
"source": [
"In the next example, we sum together the contributions for both capital growth and technology:"
]
},
{
"cell_type": "code",
"execution_count": 79,
"id": "b5747b24",
"metadata": {
"id": "b5747b24",
"outputId": "b23404ea-0275-4a57-812e-a675e5f3d37d"
},
"outputs": [
{
"data": {
"text/plain": [
"{'x2': array([[-0. , 0.0185721, 0. ],\n",
" [ 0. , 0.472 , 0. ],\n",
" [ 0. , 0. , 0.403 ]]),\n",
" 'x': array([[-0. , 0.03714422, 0.00000001],\n",
" [ 0. , 0.944 , 0. ],\n",
" [ 0. , 0. , 0.806 ]]),\n",
" 'c': array([[-3.66037866],\n",
" [-2.71090383],\n",
" [-0.07821726]]),\n",
" 'xw': array([[ 0. , 0. , 0. , 0. , 0. , 0. , 0.00399962, 0.00173897, 0. ],\n",
" [ 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.02478025, 0. ],\n",
" [ 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , -0.21392032]]),\n",
" 'w': array([[0.00799924, 0.00347793, 0. ],\n",
" [0. , 0.04956051, 0. ],\n",
" [0. , 0. , 0.42784065]]),\n",
" 'xx': array([[ 0. , -0. , 0. , -0. , -0.00021926, 0. , 0. , 0. , 0. ],\n",
" [ 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ],\n",
" [ 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.097 ]])}"
]
},
"execution_count": 79,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"lq_sum([X0_tp1, ModelSol['X1_tp1'], 0.5 * ModelSol['X2_tp1']]).coeffs"
]
},
{
"cell_type": "markdown",
"id": "5c1952e3",
"metadata": {
"id": "5c1952e3"
},
"source": [
"## 4.2 `LinQuadVar` Split and Concat \n",
"`split` breaks multiple dimensional LinQuad into one-dimensional LinQuads, while `concat` does the inverse."
]
},
{
"cell_type": "code",
"execution_count": 62,
"id": "d085149f",
"metadata": {
"id": "d085149f",
"outputId": "a6bdf64b-17a6-48c5-b159-5715a34e0174"
},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 62,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"gk_tp1, Z_tp1, Y_tp1 = ModelSol['X1_tp1'].split()\n",
"concat([gk_tp1, Z_tp1, Y_tp1])"
]
},
{
"cell_type": "markdown",
"id": "e67b04c1",
"metadata": {
"id": "e67b04c1"
},
"source": [
"## 4.3 Evaluate a `LinQuadVar` \n",
"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."
]
},
{
"cell_type": "code",
"execution_count": 72,
"id": "5782676a",
"metadata": {
"id": "5782676a",
"outputId": "0e8e4471-839b-4679-c7a2-fd8d6a8e6879"
},
"outputs": [
{
"data": {
"text/plain": [
"array([[ -3.66764679],\n",
" [ -2.70809391],\n",
" [ 0.00945342],\n",
" [ -0.03677232],\n",
" [-12.41358991]])"
]
},
"execution_count": 72,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x1 = np.zeros([n_X ,1])\n",
"x2 = np.zeros([n_X ,1])\n",
"w = np.random.multivariate_normal(np.zeros(n_W),np.eye(n_W),size = 1).T\n",
"ModelSol['JX_tp1'](*(x1,x2,w))"
]
},
{
"cell_type": "markdown",
"id": "5051f96c",
"metadata": {
"id": "5051f96c"
},
"source": [
"## 4.4 Next period expression for `LinQuadVar` \n",
"`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:"
]
},
{
"cell_type": "code",
"execution_count": 73,
"id": "7e088a47",
"metadata": {
"id": "7e088a47",
"outputId": "116af9b7-9791-4391-b7a4-8d9ee133baff"
},
"outputs": [
{
"data": {
"text/latex": [
"$\\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}$"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"cmk1_tp1 = next_period(ModelSol['J1_t'][0], ModelSol['X1_tp1'])\n",
"disp(cmk1_tp1, '\\\\log\\\\frac{C_{t+\\epsilon}^1}{K_{t+\\epsilon}^1}') "
]
},
{
"cell_type": "code",
"execution_count": 74,
"id": "fa8ab664",
"metadata": {
"id": "fa8ab664",
"outputId": "f4244fa2-2a8a-4a42-8a33-2ab7407c7c8d"
},
"outputs": [
{
"data": {
"text/latex": [
"$\\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}$"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"cmk2_tp1 = next_period(ModelSol['J2_t'][0], ModelSol['X1_tp1'], ModelSol['X2_tp1'])\n",
"disp(cmk2_tp1, '\\\\log\\\\frac{C_{t+\\epsilon}^2}{K_{t+\\epsilon}^2}') "
]
},
{
"cell_type": "markdown",
"id": "daa0c4c1",
"metadata": {
"id": "daa0c4c1"
},
"source": [
"## 4.6 Compute the Expectation of time $t+\\epsilon$ `LinQuadVar` \n",
"\n",
"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:\n"
]
},
{
"cell_type": "code",
"execution_count": 75,
"id": "3140272c",
"metadata": {
"id": "3140272c",
"outputId": "418a75e2-ff56-4bd7-8471-67634f26654c"
},
"outputs": [
{
"data": {
"text/latex": [
"$\\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}$"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"E_w = ModelSol['util_sol']['μ_0']\n",
"cov_w = np.eye(n_W)\n",
"E_ww = cal_E_ww(E_w, cov_w)\n",
"E_cmk2_tp1 = E(cmk2_tp1, E_w, E_ww)\n",
"disp(E_cmk2_tp1, '\\mathbb{E}[\\\\log\\\\frac{C_{t+\\epsilon}^2}{K_{t+\\epsilon}^2}|\\mathfrak{F_t}]')"
]
},
{
"cell_type": "markdown",
"id": "18a1686c",
"metadata": {
"id": "18a1686c"
},
"source": [
"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`."
]
},
{
"cell_type": "code",
"execution_count": 76,
"id": "a58482f7",
"metadata": {
"id": "a58482f7",
"outputId": "b8f2fc85-7c37-4d01-a5ed-c6b4d6d940a0"
},
"outputs": [
{
"data": {
"text/latex": [
"$\\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}$"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"N_cm = N_tilde_measure(ModelSol['util_sol']['log_N_tilde'],(1,n_X,n_W))\n",
"E_cmk2_tp1_tilde = E_N_tp1(cmk2_tp1, N_cm)\n",
"disp(E_cmk2_tp1_tilde, '\\mathbb{\\\\tilde{E}}[\\\\log\\\\frac{C_{t+\\epsilon}^2}{K_{t+\\epsilon}^2}|\\mathfrak{F_t}]')"
]
}
],
"metadata": {
"colab": {
"provenance": []
},
"interpreter": {
"hash": "1c67a9c81f4b1f676ee52b818d3e176be5c9dd44f82a35d6054954eb94c37bc5"
},
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.5"
}
},
"nbformat": 4,
"nbformat_minor": 5
}