{
"cells": [
{
"cell_type": "markdown",
"id": "a985108b-129f-4f61-b092-a3392e4a0daf",
"metadata": {},
"source": [
"# 6. Heterogeneous Agent Model in Continous Time"
]
},
{
"cell_type": "markdown",
"id": "b11fe52c-6710-4ad2-aa10-c0af2c0979fe",
"metadata": {},
"source": [
"DING Minjie, Spring 2025"
]
},
{
"cell_type": "markdown",
"id": "2e407648-7067-412e-a1fc-28146c96d3f9",
"metadata": {},
"source": [
"__________________________________________________________________________________________________________________________________________"
]
},
{
"cell_type": "markdown",
"id": "22b6bc49-3039-4829-b915-17b91d1ed805",
"metadata": {},
"source": [
"We will discuss Aiyagari model in continous time in this tutorial. The main reference is \"*Income and Wealth Distribution in Macroeconomics: A Continuous-Time Approach*\".\n",
"\n",
"1. Basic Idea in Continous Time\n",
"2. Income in Ornstein-Uhlenbeck Process\n",
"3. Evolution Matrix of Productivity (a Ornstein-Uhlenbeck Process)\n",
"4. Evolution Matrix of Asset\n",
"5. Hamilton–Jacobi–Bellman equation\n",
"6. Kolmogorov Forward Equation/ Fokker-Planck Equation\n",
"7. Income in Poisson Process"
]
},
{
"cell_type": "markdown",
"id": "89a1b398-a3dd-48e3-8809-833f28bb3359",
"metadata": {},
"source": [
"________________________________________________________________________________________________________________________________________________"
]
},
{
"cell_type": "markdown",
"id": "ad5c837a-4829-4ff1-8001-2a03b995d3ae",
"metadata": {},
"source": [
"# 1. Basic Idea in Continous Time"
]
},
{
"cell_type": "markdown",
"id": "9b0d970e-ee84-42ac-b836-2ae097fabaf4",
"metadata": {},
"source": [
"_______________________________________________________________________________________________________________________________________________"
]
},
{
"cell_type": "markdown",
"id": "d30f4338-5417-4900-ab4a-6b86bd54a5d4",
"metadata": {},
"source": [
"In discrete case, we have income grid and income transition on one dimension. On the other dimension, we have asset grid, combined with income process, we get: asset transition (asset policy function), and know how asset and consumption affect value function.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "140218a0-5a43-40e8-a914-36b989781d14",
"metadata": {},
"source": [
"Similarly, in continous case, we need both. \n",
"\n",
"On the one dimension, we need income grid and its transition. Sometimes, income is determined by labor productivity, so we use productivity $z$ here, in consistent with literature. Two popular methods to construct continous income are: Poisson Process and Ornstein-Uhlenbeck Process. \n",
"\n",
"On the other dimension, we need to know how productivity and asset affect value function. In continous case, we use Differential. \n",
" \n",
"When asset increase from low to high, value function is expected to increase, and we use forward difference to approximate it. That is, use high value function and current value function difference to approximate the increasing process:\n",
" $$ \n",
" v'_j(a_i) \\approx \\frac{v_{i+1,j} - v_{i,j}}{\\Delta a}\n",
" $$\n",
"\n",
"When asset decrease from high to low, value function is expected to decrease, we use backward difference to approximate it. That is, use current value and a low value to approximate the decreasing process:\n",
" $$\n",
" v'_j(a_i) \\approx \\frac{v_{i,j} - v_{i-1,j}}{\\Delta a}\n",
" $$"
]
},
{
"cell_type": "markdown",
"id": "34e4b62b-ab87-4f38-8f65-f105132a846a",
"metadata": {},
"source": [
"It the same for productivity.\n",
" \n",
"When productivity increase from low to high, value function is expected to increase, and we use forward difference to approximate it:\n",
" $$ \n",
" v'_i(z_j) \\approx \\frac{v_{i,j+1} - v_{i,j}}{\\Delta z}\n",
" $$\n",
"\n",
"When productivity decrease from high to low, value function is expected to decrease, we use backward difference to approximate it:\n",
" $$\n",
" v'_i(z_j) \\approx \\frac{v_{i,j} - v_{i,j-1}}{\\Delta z}\n",
" $$"
]
},
{
"cell_type": "markdown",
"id": "0b6d7d84-c494-4c7b-b1d5-21771c1a001d",
"metadata": {},
"source": [
"Note the subscript here, $i$ is for different asset $a$, $j$ is for different productivity $z$. "
]
},
{
"cell_type": "markdown",
"id": "67536af8-a5df-4339-977e-ed1d9e78ee72",
"metadata": {},
"source": [
"_________________________________________________________________________________________________________________________________________________"
]
},
{
"cell_type": "markdown",
"id": "3bda1ac9-2d0a-4463-b0b5-7a23c56cd0e0",
"metadata": {},
"source": [
"# 2. Income in Ornstein-Uhlenbeck Process"
]
},
{
"cell_type": "raw",
"id": "f28709c8-6240-4115-8007-fa8f5b7bc791",
"metadata": {},
"source": [
"___________________________________________________________________________________________________________________________________________________"
]
},
{
"cell_type": "markdown",
"id": "b1747335-cf1e-4dd9-8004-8fb253b7aaaa",
"metadata": {},
"source": [
"The Ornstein-Uhlenbeck (O-U) process is a type of stochastic process often used to model mean-reverting behavior. It can be described by the following stochastic differential equation (SDE):\n",
"\n",
"$$\n",
"dX_t = \\theta(\\mathbb{E}(X_t) - X_t)dt + \\sigma dW_t\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "3cc75705-d291-4f69-9f84-81785616cb5e",
"metadata": {},
"source": [
"where:\n",
"- $X_t$ is the value of the process at time $t$.\n",
"- $dW_t$ represents the increment of a Wiener process (or Brownian motion).\n"
]
},
{
"cell_type": "markdown",
"id": "afc4d8b8-5b9e-4707-b9bf-a30af84b4f98",
"metadata": {},
"source": [
"- $E(X_t)$: The long-term mean level to which the process reverts.\n",
"- $\\theta$: The speed of mean reversion, indicating how quickly the process returns to the mean after a disturbance.\n",
"- $\\sigma$: The volatility of the process, representing the intensity of random fluctuations.\n"
]
},
{
"cell_type": "markdown",
"id": "d43af476-1f6b-41a1-8a31-3104e178e3a0",
"metadata": {},
"source": [
"The O-U process can be solved analytically. The solution to the SDE is given by:\n",
"\n",
" $$ \n",
" X_t = X_0 e^{-\\theta t} + \\mathbb{E}(X_t)(1 - e^{-\\theta t}) + \\sigma \\int_0^t e^{-\\theta(t-s)} dW_s \n",
" $$\n",
"\n",
"Here, $X_0$ is the initial value of the process.\n"
]
},
{
"cell_type": "markdown",
"id": "d8848c6d-d3ab-4abf-a346-0d44a9a4e7f7",
"metadata": {},
"source": [
"To find the variance, we need to calculate $\\text{Var}(X_t)$, which is defined as:\n",
"\n",
" $$\n",
" \\text{Var}(X_t) = \\mathbb{E}[(X_t - \\mathbb{E}[X_t])^2] \n",
" $$\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "a8ef82b8-ba20-4637-8098-a6bc0bbbe7b7",
"metadata": {},
"source": [
"The stochastic integral term $\\int_0^t e^{-\\theta(t-s)} dW_s$ is a key component. The variance of this term can be calculated using the properties of Itô integrals:\n",
"\n",
" $$\n",
" \\text{Var}\\left(\\int_0^t e^{-\\theta(t-s)} dW_s\\right) = \\int_0^t e^{-2\\theta(t-s)} ds \n",
" $$\n",
"\n",
" Evaluating this integral gives:\n",
"\n",
" $$ \n",
" = \\frac{1}{2\\theta}(1 - e^{-2\\theta t}) \n",
" $$\n"
]
},
{
"cell_type": "markdown",
"id": "839a9165-e459-47b6-98a0-d76c2258cf83",
"metadata": {},
"source": [
"Hence, the variance of $X_t$ can be derived as:\n",
"\n",
" $$\n",
" \\text{Var}(X_t) = \\sigma^2 \\frac{1 - e^{-2\\theta t}}{2\\theta} \n",
" $$\n",
"\n",
" As $t \\to \\infty$, the term $e^{-2\\theta t}$ approaches zero, leading to the steady-state variance:\n",
"\n",
" $$ \n",
" \\text{Var}(X_t) = \\frac{\\sigma^2}{2\\theta} \n",
" $$\n"
]
},
{
"cell_type": "markdown",
"id": "d7a16f7b-8315-4a20-9574-9a23f99763b4",
"metadata": {},
"source": [
"__________________________________________________________________________________________________________________________________________"
]
},
{
"cell_type": "markdown",
"id": "24fae498-2820-4f14-93dc-fd82dec8e75b",
"metadata": {},
"source": [
"# 3. Evolution Matrix of Productivity (a Ornstein-Uhlenbeck Process)"
]
},
{
"cell_type": "markdown",
"id": "0afd4d7e-bb59-49c9-a2b2-2f607f32a66c",
"metadata": {},
"source": [
"____________________________________________________________________________________________________________________________________________"
]
},
{
"cell_type": "markdown",
"id": "70c54aea-6076-46d0-a506-41b2129c3bb8",
"metadata": {},
"source": [
"Now suppose the productivity $z_t$ follows Ornstein-Uhlenbeck process.\n",
"\n",
"$$\n",
"dz_t = \\theta (\\mathbb{E}(z_t) - z_t)dt + \\sigma dW_t\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "476342ce-463d-4b5d-a131-1a028d8c36fc",
"metadata": {},
"source": [
"For the convenience of expression, define the drift term as:\n",
"\n",
"$$\n",
"\\mu = \\theta (\\mathbb{E}(z_t)-z_t)\n",
"$$\n",
"\n",
"Then the Ornstein-Uhlenbeck process can be written as:\n",
"\n",
"$$\n",
"dz_t = \\mu dt + \\sigma dW_t\n",
"$$\n",
"\n",
"For the convenience of expression, we call variance as:\n",
"\n",
"$$\n",
"s^2 = \\sigma^2\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "f6561bd8-6ede-4592-8ec7-a0f76d3673f1",
"metadata": {},
"source": [
"**1. When $z_t$ is larger than $\\mathbb{E}(z_t)$**, \n",
"\n",
"$\\mu = \\theta (\\mathbb{E}(z_t)-z_t) < 0$, \n",
"\n",
"$dz_t = \\mu dt + \\sigma dW_t$ has high possibility to be negative, and $z_t$ tends to decrease from high to low."
]
},
{
"cell_type": "markdown",
"id": "4ed45dac-9ad2-4dc5-b30b-98d407c9af1c",
"metadata": {},
"source": [
"Define **Chi ($\\chi $)**:\n",
"\n",
"\n",
"$$\n",
"\\chi = -\\frac{\\min(\\mu, 0)}{\\Delta z} + \\frac{s^2}{2 \\Delta z^2}\n",
"$$\n",
"\n",
"The first term, $-\\frac{\\min(\\mu, 0)}{\\Delta z}$, when productivity $ z $ is above its mean, this term is positively large.\n",
"\n",
"When put it in lower diagonal, it reflect tendency to change from higher to lower index of grid, or from high productivity to lower productivity, that is to revert back to the mean. "
]
},
{
"cell_type": "markdown",
"id": "ef3f9ab7-0180-40ac-8017-6f4f72a6ee91",
"metadata": {},
"source": [
"**2. When $z_t$ is smaller than $\\mathbb{E}(z_t)$**, \n",
"\n",
"$\\mu = \\theta (\\mathbb{E}(z_t)-z_t) > 0$, \n",
"\n",
"$dz_t = \\mu dt + \\sigma dW_t$ has high possibility to be positive, and $z_t$ tends to increase from low to high."
]
},
{
"cell_type": "markdown",
"id": "c5d8cca2-fc89-4089-89bd-4fdb9046faaa",
"metadata": {},
"source": [
"Define **Zeta ($ \\zeta $)**:\n",
"\n",
"$$\n",
"\\zeta = \\frac{\\max(\\mu, 0)}{\\Delta z} + \\frac{s^2}{2 \\Delta z^2}\n",
"$$\n",
"\n",
"When put upper lower diagonal, it reflect tendency to change from low to high index of grid, or from lower productivity to higher productivity, that is to revert back to the mean. "
]
},
{
"cell_type": "markdown",
"id": "dc137266-b4db-496d-9965-8c5a79fc592a",
"metadata": {},
"source": [
"**3. When $z_t$ is close to $\\mathbb{E}(z_t) $**:\n",
"\n",
"Suppose the difference interval is not small enough to capture the difference between $z_t$ and $\\mathbb{E}(z_t)$, $z_t$ stay at the original state."
]
},
{
"cell_type": "markdown",
"id": "44f6dd8d-765f-46ff-947d-5d64fbf6e034",
"metadata": {},
"source": [
"Define **YY ($ \\text{yy} $)**\n",
"\n",
"$$\n",
"\\text{yy} = -\\chi - \\zeta = \\frac{\\min(\\mu, 0)}{\\Delta z} - \\frac{\\max(\\mu, 0)}{\\Delta z} - \\frac{s^2}{\\Delta z^2}\n",
"$$\n",
"\n",
"When put it in the diagonal, it means $z_t$ stay at the original state."
]
},
{
"cell_type": "markdown",
"id": "9ae04628-080d-4b0b-8515-6a2ddb71d68e",
"metadata": {},
"source": [
"Finally, we can construct a sparse matrix to show transition of productivity. \n",
"\n",
"1. **Chi ($ \\chi $)**: used for the lower diagonal of the sparse matrix.\n",
"\n",
"2. **YY ($ \\text{yy} $)**: forms the main diagonal of the sparse matrix.\n",
"\n",
"3. **Zeta ($ \\zeta $)**: used for the upper diagonal of the sparse matrix."
]
},
{
"cell_type": "markdown",
"id": "c06f8a17-0e0c-4461-86c8-ec996a8ca519",
"metadata": {},
"source": [
"Remember the **spdiags** function introduced before, we can use 'spdiags' to get sparse matrix."
]
},
{
"cell_type": "markdown",
"id": "85540c6a-c4b2-40a8-84d4-882dfb8199bd",
"metadata": {},
"source": [
"Before construct the sparse matrix, we set parameters, and the grids. The following three code cells are:\n",
" \n",
" Parameters\n",
" Grids\n",
" Sparse Matrix for Productivity"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "754f15ac-3531-493b-ba12-32d0418dba4b",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"ans = 'C:\\Users\\ading\\A_TA_Notes'
"
],
"text/plain": [
"ans = 'C:\\Users\\ading\\A_TA_Notes'"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"cd\"C:\\Users\\ading\\A_TA_Notes\"\n",
"pwd"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "76caa371-2cc3-480a-b7eb-5c63cb0ce528",
"metadata": {},
"outputs": [],
"source": [
"clear all; \n",
"close all; \n",
"format long; % Set long format for outputs\n",
"tic; % Start timing the execution"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "e1259bad-70b0-47bc-8b16-7ad1640cc372",
"metadata": {},
"outputs": [],
"source": [
"%% PARAMETERS\n",
"% household and firm\n",
"ga = 2; % CRRA utility with parameter gamma\n",
"rho = 0.05; % Discount rate\n",
"alpha = 0.35; % Capital share of production function F = K^alpha * L^(1-alpha) \n",
"delta = 0.1; % Capital depreciation rate\n",
"\n",
"% productivity and O-U process\n",
"zmean = 1.0; % Mean of the Ornstein-Uhlenbeck (O-U) process (in levels)\n",
"sig2 = (0.10)^2; % Variance of the O-U process (sigma^2)\n",
"Corr = exp(-0.3); % Persistence of the O-U process (-log(Corr))\n",
"the = -log(Corr); % Calculate the rate of mean reversion\n",
"Var = sig2/(2*the); % Variance for the O-U process\n",
"\n",
"\n",
"% Iteration parameters\n",
"maxit = 100; % Maximum number of iterations in the HJB loop\n",
"maxitK = 100; % Maximum number of iterations in the K loop\n",
"crit = 10^(-6); % Convergence criterion for HJB loop\n",
"critK = 1e-5; % Convergence criterion for K loop\n",
"Delta = 1000; % Delta in HJB algorithm\n",
"K = 3.8; % Initial aggregate capital; should be close to the solution for convergence\n",
"relax = 0.99; % Relaxation parameter for the algorithm"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "204c7a00-50e7-41c4-8c3f-94473639b22d",
"metadata": {},
"outputs": [],
"source": [
"%% Grid\n",
"% grid dimension\n",
"J = 40; % Number of points for productivity (z)\n",
"I = 100; % Number of points for wealth (a)\n",
"\n",
"% productivity \n",
"zmin = 0.5; % Minimum value of z\n",
"zmax = 1.5; % Maximum value of z\n",
"z = linspace(zmin, zmax, J); % Productivity vector (row vector of productivity values)\n",
"dz = (zmax - zmin) / (J - 1); % Step size for productivity\n",
"dz2 = dz^2; % Square of the step size for productivity\n",
"zz = ones(I, 1) * z; % Replicate productivity vector for matrix operations\n",
"\n",
"% asset\n",
"amin = -1; % Borrowing constraint (minimum wealth)\n",
"amax = 30; % Maximum wealth\n",
"a = linspace(amin, amax, I)'; % Wealth vector (column vector of wealth values)\n",
"da = (amax - amin) / (I - 1); % Step size for wealth\n",
"aa = a * ones(1, J); % Replicate wealth vector for matrix operations\n",
"\n",
"% Preallocation\n",
"% Finite difference approximation of the partial derivatives\n",
"Vaf = zeros(I, J); % Forward difference value function\n",
"Vab = zeros(I, J); % Backward difference value function\n",
"Vzf = zeros(I, J); % Forward difference for z\n",
"Vzb = zeros(I, J); % Backward difference for z\n",
"Vzz = zeros(I, J); % Second derivative for z\n",
"c = zeros(I, J); % Consumption matrix"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "82c385a3-7951-4850-b7f6-3280e3a8328a",
"metadata": {},
"outputs": [],
"source": [
"%% Construct Matrix Aswitch summarizing evolution of z\n",
"mu = the * (zmean - z); % Drift term (calculated from Itô's lemma)\n",
"s2 = sig2 .* ones(1, J); % Variance term (constant for all z)\n",
"\n",
"chi = -min(mu, 0) / dz + s2 / (2 * dz2); % Coefficient for upward diffusion\n",
"yy = min(mu, 0) / dz - max(mu, 0) / dz - s2 / dz2; % Coefficient for the drift term\n",
"zeta = max(mu, 0) / dz + s2 / (2 * dz2); % Coefficient for downward diffusion\n",
"\n",
"% This will be the upper diagonal of the matrix Aswitch\n",
"updiag = zeros(I, 1); % Initialize upper diagonal\n",
"for j = 1:J\n",
" updiag = [updiag; repmat(zeta(j), I, 1)]; % Fill upper diagonal with zeta values\n",
"end\n",
"\n",
"% This will be the center diagonal of the matrix Aswitch\n",
"centdiag = repmat(chi(1) + yy(1), I, 1); % Initialize center diagonal\n",
"for j = 2:J-1\n",
" centdiag = [centdiag; repmat(yy(j), I, 1)]; % Fill center diagonal with yy values\n",
"end\n",
"centdiag = [centdiag; repmat(yy(J) + zeta(J), I, 1)]; % Add last entry\n",
"\n",
"% This will be the lower diagonal of the matrix Aswitch\n",
"lowdiag = repmat(chi(2), I, 1); % Initialize lower diagonal\n",
"for j = 3:J\n",
" lowdiag = [lowdiag; repmat(chi(j), I, 1)]; % Fill lower diagonal with chi values\n",
"end\n",
"\n",
"% Add up the upper, center, and lower diagonal into a sparse matrix\n",
"Aswitch = spdiags(centdiag, 0, I*J, I*J) + spdiags(lowdiag, -I, I*J, I*J) + spdiags(updiag, I, I*J, I*J);"
]
},
{
"cell_type": "markdown",
"id": "123ed9d5-5b41-4277-b462-39eecfddbba1",
"metadata": {},
"source": [
"_______________________________________________________________________________________________________________________________________________"
]
},
{
"cell_type": "markdown",
"id": "388e459f-209c-4bd8-a1b5-354c49504e32",
"metadata": {},
"source": [
"# 4. Evolution Matrix of Asset"
]
},
{
"cell_type": "markdown",
"id": "b6ead0a9-ce1d-4d93-8693-5a18b44b64ac",
"metadata": {},
"source": [
"__________________________________________________________________________________________________________________________________________________"
]
},
{
"cell_type": "markdown",
"id": "8e6690d9-d52b-45dd-8518-f4e205ddbcc7",
"metadata": {},
"source": [
"Given a guessed value function, we can construct evolution matrix of asset, similar to that of productivity."
]
},
{
"cell_type": "markdown",
"id": "bb8a5ba7-6349-4837-b720-134bfaa58119",
"metadata": {},
"source": [
"An important distinction between continous case and discrete case lie in FOC condition."
]
},
{
"cell_type": "markdown",
"id": "574de831-b08d-4ada-9e0e-70e777fd0163",
"metadata": {},
"source": [
"Assume CRRA utility, so the two first order conditions are\n",
"\n",
"**continuous time**\n",
"\n",
"$$\n",
"c^{-\\gamma} = v'_j(a), j = 1,2 \n",
"$$\n",
"\n",
"**discrete time**\n",
"\n",
"$$\n",
"c^{-\\gamma} ≥ \\beta \\sum_{k=1}^{2} \\pi_{jk} v'_k(a'), \\quad a' = y_j +(1+r)a-c,\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "02cd8e60-adca-4668-932d-700c45b3556c",
"metadata": {},
"source": [
"The first-order condition of continous time is “static” in the sense that it only involves contemporaneous variables. Given (a guess for) the value function $v_j(a)$ it can be solved by hand. In contrast, the discrete-time condition defines the optimal choice only implicitly."
]
},
{
"cell_type": "markdown",
"id": "cd123514-6925-48df-94d3-3de07b5564e8",
"metadata": {},
"source": [
"From the FOC condition, we can get $c$, since value function is given.\n",
"\n",
"Then we can get saving $s$ from budget constraint.\n",
"\n",
"And $s$ is exactly the change of asset, or \"drift item\" of asset, if analogous to productivity.\n",
"\n",
"Then, construct the transition matrix similar to that of productivity."
]
},
{
"cell_type": "markdown",
"id": "8eeb4abd-0eb3-4493-82ea-bf7d8e85b043",
"metadata": {},
"source": [
"__________________________________________________________________________________________________________________________________________"
]
},
{
"cell_type": "markdown",
"id": "cfc72831-1969-444e-bbbd-11d4bb0385ab",
"metadata": {},
"source": [
"# 5. Hamilton–Jacobi–Bellman equation"
]
},
{
"cell_type": "markdown",
"id": "455b661d-f993-4ecf-82c3-c161e2e2ff81",
"metadata": {},
"source": [
"_____________________________________________________________________________________________________________________________________________"
]
},
{
"cell_type": "markdown",
"id": "9785f4be-c433-4c62-996e-8b1a4fa03014",
"metadata": {},
"source": [
"The HJB equation in this case is:"
]
},
{
"cell_type": "markdown",
"id": "43fee2a4-205f-4dbe-adc3-3d5c0cd953c9",
"metadata": {},
"source": [
"$$\n",
"\\frac{v_{i, j}^{n+1}-v_{i, j}^n}{\\Delta}+\\rho v_{i, j}^{n+1}= u\\left(c_{i, j}^n\\right)+\\frac{v_{i+1, j}^{n+1}-v_{i, j}^{n+1}}{\\Delta a}\\left(s_{i, j, F}^n\\right)^{+}+\\frac{v_{i, j}^{n+1}-v_{i-1, j}^{n+1}}{\\Delta a}\\left(s_{i, j, B}^n\\right)^{-} +\\frac{v_{i, j+1}^{n+1}-v_{i, j}^{n+1}}{\\Delta z} \\mu_j^{+}+\\frac{v_{i, j}^{n+1}-v_{i, j-1}^{n+1}}{\\Delta z} \\mu_j^{-}+\\frac{\\sigma_j^2}{2} \\frac{v_{i, j+1}^{n+1}-2 v_{i, j}^{n+1}+v_{i, j-1}^{n+1}}{(\\Delta z)^2}\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "7b4546aa-47aa-48fc-ba86-95e953567c61",
"metadata": {},
"source": [
"$\\frac{1}{\\Delta}$ : a very small interval "
]
},
{
"cell_type": "markdown",
"id": "d6cb441b-bd44-422f-96bf-77f8cd4f9b92",
"metadata": {},
"source": [
"$v^{n}_{i,j}$ : old guessed value function"
]
},
{
"cell_type": "markdown",
"id": "e09ba6a8-8c97-45a9-91bf-8aeff22a26b8",
"metadata": {},
"source": [
"$v^{n+1}_{i,j}$: new value function from iteration"
]
},
{
"cell_type": "markdown",
"id": "2cb407fd-dba3-44ff-a464-be830bf46c85",
"metadata": {},
"source": [
"$\\rho$ : discount rate"
]
},
{
"cell_type": "markdown",
"id": "14b13e8b-1f00-4730-b055-f3304993bb64",
"metadata": {},
"source": [
"$\\frac{v_{i+1, j}^{n+1}-v_{i, j}^{n+1}}{\\Delta a}$ : forward difference with respect to asset"
]
},
{
"cell_type": "markdown",
"id": "b6848684-5fe1-4739-bf8a-4f2290af5e5f",
"metadata": {},
"source": [
"$(s_{i, j, F}^n)^{+}$ : positive change of asset (positive saving)"
]
},
{
"cell_type": "markdown",
"id": "99142634-651a-4233-b875-2ae47d7e4450",
"metadata": {},
"source": [
"$\\frac{v_{i+1, j}^{n+1}-v_{i, j}^{n+1}}{\\Delta a}(s_{i, j, F}^n)^{+}$ : value gained from increasing saving"
]
},
{
"cell_type": "markdown",
"id": "d313141b-0485-4072-aa94-663e6ef9b006",
"metadata": {},
"source": [
"$\\frac{v_{i, j}^{n+1}-v_{i-1, j}^{n+1}}{\\Delta a}(s_{i, j, B}^n)^{-}$ : value gained from decreasing saving"
]
},
{
"cell_type": "markdown",
"id": "13df7d9e-b856-4e83-8225-1fd6b27610d8",
"metadata": {},
"source": [
"$\\frac{v_{i+1, j}^{n+1}-v_{i, j}^{n+1}}{\\Delta z}(\\mu_{i, j, F}^n)^{+}$ : value gained from increasing productivity"
]
},
{
"cell_type": "markdown",
"id": "53e055e4-9894-4f67-adbc-61340017b303",
"metadata": {},
"source": [
"$\\frac{v_{i, j}^{n+1}-v_{i-1, j}^{n+1}}{\\Delta z}(\\mu_{i, j, B}^n)^{-}$ : value gained from decreasing productivity"
]
},
{
"cell_type": "markdown",
"id": "9c2469f8-54e4-4fde-88cc-ae543b3c9f9e",
"metadata": {},
"source": [
"$\\frac{\\sigma_j^2}{2} \\frac{v_{i, j+1}^{n+1}-2 v_{i, j}^{n+1}+v_{i, j-1}^{n+1}}{(\\Delta z)^2}$ : value gained from volatility"
]
},
{
"cell_type": "markdown",
"id": "9be8b8af-264f-478a-94b2-de2b9b84fef1",
"metadata": {},
"source": [
"Rearrange HJB equation, and collect terms with the same subscripts on the right-hand side."
]
},
{
"cell_type": "markdown",
"id": "7f280812-c948-4d08-a712-287283c1efbf",
"metadata": {},
"source": [
"$$\n",
"\\begin{aligned}\n",
"& \\frac{v_{i, j}^{n+1}-v_{i, j}^n}{\\Delta}+\\rho v_{i, j}^{n+1}=u\\left(c_{i, j}^n\\right)+v_{i-1, j}^{n+1} x_{i, j}+v_{i, j}^{n+1}\\left(y_{i, j}+v_j\\right)+v_{i+1, j}^{n+1} z_{i, j}+v_{i, j-1}^{n+1} \\chi_j+v_{i, j+1}^{n+1} \\zeta_j \\\\\n",
"& x_{i, j}=-\\frac{\\left(s_{i, j, B}^n\\right)^{-}}{\\Delta a} \\\\\n",
"& y_{i, j}=-\\frac{\\left(s_{i, j, F}^n\\right)^{+}}{\\Delta a}+\\frac{\\left(s_{i, j, B}^n\\right)^{-}}{\\Delta a} \\\\\n",
"& z_{i, j}=\\frac{\\left(s_{i, j, F}^n\\right)^{+}}{\\Delta a} \\\\\n",
"& \\chi_j=-\\frac{\\mu_j^{-}}{\\Delta z}+\\frac{\\sigma_j^2}{2(\\Delta z)^2} \\\\\n",
"& v_j=\\frac{\\mu_j^{-}}{\\Delta z}-\\frac{\\mu_j^{+}}{\\Delta z}-\\frac{\\sigma_j^2}{(\\Delta z)^2} \\\\\n",
"& \\zeta_j=\\frac{\\mu_j^{+}}{\\Delta z}+\\frac{\\sigma_j^2}{2(\\Delta z)^2}\n",
"\\end{aligned}\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "89b9084b-91df-43c8-aaa2-dcb262eeea72",
"metadata": {},
"source": [
"Let:\n",
"$$\n",
"AA = x_{i,j} + y_{i,j} + z_{i,j}\n",
"$$\n",
"$$\n",
"Aswitch = \\chi_{j} + v_j + \\zeta_j\n",
"$$\n",
"$$\n",
"A = AA + Aswitch\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "28c95a82-4650-4adc-b9dc-7106ee07eda5",
"metadata": {},
"source": [
"Then put $v^n_{i,j}$ and $v^{n+1}_{i,j}$ on different hand, we get: "
]
},
{
"cell_type": "markdown",
"id": "dd8612a0-3c09-472a-bf41-542c35279dbb",
"metadata": {},
"source": [
"$$\n",
"(\\frac{1}{\\Delta} + \\rho - A) v^{n+1}_{i,j} = u(c^n_{i,j}) + \\frac{1}{\\Delta} v^n_{i,j}\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "6d355784-dbe5-4850-8891-3c9c8a4106e1",
"metadata": {},
"source": [
"Given guessed value function $v^n_{i,j}$ and solved $c_{i,j}^n$ from it, we can then get new value function.\n",
"\n",
"Then iterate, until value function converge. "
]
},
{
"cell_type": "markdown",
"id": "b7a305cf-e6b7-4b08-a068-2375e131796c",
"metadata": {},
"source": [
"Let:\n",
"$$\n",
"b = u(c^n_{i,j}) + \\frac{1}{\\Delta} v^n_{i,j}\n",
"$$\n",
"$$\n",
"B = (\\frac{1}{\\Delta} + \\rho - A)\n",
"$$\n",
"Then: \n",
"$$\n",
"B v^{n+1}_{i,j} = b\n",
"$$\n",
"$$\n",
"v^{n+1}_{i,j} = B^{-1} b\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "f5f77e02-7d6c-498e-b363-23d6b8f78a6e",
"metadata": {},
"source": [
"Until:"
]
},
{
"cell_type": "markdown",
"id": "6b67cc29-924e-48f4-bf22-00b36db9e574",
"metadata": {},
"source": [
"$$\n",
"abs (v^{n+1}_{i,j} - v^{n}_{i,j}) \\to \\epsilon\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "0ab91ae1-7414-4e99-878a-73765155cf6f",
"metadata": {},
"source": [
"___________________________________________________________________________________________________________________________________________"
]
},
{
"cell_type": "markdown",
"id": "839dfaae-816a-4e0e-8e4c-ca195231ded8",
"metadata": {},
"source": [
"# 6. Kolmogorov Forward Equation/ Fokker-Planck Equation"
]
},
{
"cell_type": "markdown",
"id": "97087b05-966e-4e88-bb64-53f055a3f29e",
"metadata": {},
"source": [
"____________________________________________________________________________________________________________________________________________"
]
},
{
"cell_type": "markdown",
"id": "798310a8-475f-4681-9370-4901838ae286",
"metadata": {},
"source": [
"The Kolmogorov Equation describes the time evolution of the probability density function of a stochastic process."
]
},
{
"cell_type": "markdown",
"id": "e8967e28-1e56-4dd8-8d16-57d44e3600f4",
"metadata": {},
"source": [
"In stationary equilibrium, changes from productivity or saving will not affect the distribution of asset, (will not change probability density function of asset). Let $gg$ be the distribution, then: "
]
},
{
"cell_type": "markdown",
"id": "938b8855-a6b9-47fd-a6b4-a5769cd5efbe",
"metadata": {},
"source": [
"$$\n",
"gg * A' = 0\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "42ee28f9-3ef9-43be-89e5-22a18e33f6ce",
"metadata": {},
"source": [
"where $A$ is what we defined above, denoting effect on value function from productivity change and saving change. "
]
},
{
"cell_type": "markdown",
"id": "72b9d4ff-e1b4-47b8-aee5-25135762f1cd",
"metadata": {},
"source": [
"Here comes the main loop code."
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "b56cc34e-cd79-418d-ae66-ec0d7a1202f7",
"metadata": {},
"outputs": [],
"source": [
"%INITIAL GUESS\n",
"r = alpha * K^(alpha-1) -delta; %interest rates\n",
"w = (1-alpha) * K^(alpha); %wages\n",
"v0 = (w*zz + r.*aa).^(1-ga)/(1-ga)/rho;\n",
"v = v0;\n",
"dist = zeros(1,maxit);"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "13956079-fc30-4060-8a42-a27024f34629",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Main loop iteration\n",
" 1\n",
"\n",
"Value Function Converged, Iteration = \n",
" 6\n",
"\n",
"历时 132.078659 秒。\n",
" 1.460525943026954\n",
"\n",
"Main loop iteration\n",
" 2\n",
"\n",
"Value Function Converged, Iteration = \n",
" 4\n",
"\n",
"历时 132.195904 秒。\n",
" 2.020815936779077\n",
"\n",
"Main loop iteration\n",
" 3\n",
"\n",
"Value Function Converged, Iteration = \n",
" 4\n",
"\n",
"历时 132.340333 秒。\n",
" 2.660255904062295\n",
"\n",
"Main loop iteration\n",
" 4\n",
"\n",
"Value Function Converged, Iteration = \n",
" 4\n",
"\n",
"历时 132.462327 秒。\n",
" 3.225918298182704\n",
"\n",
"Main loop iteration\n",
" 5\n",
"\n",
"Value Function Converged, Iteration = \n",
" 4\n",
"\n",
"历时 132.586729 秒。\n",
" 3.561827701523755\n",
"\n",
"Main loop iteration\n",
" 6\n",
"\n",
"Value Function Converged, Iteration = \n",
" 4\n",
"\n",
"历时 132.726711 秒。\n",
" 3.691463823642045\n",
"\n",
"Main loop iteration\n",
" 7\n",
"\n",
"Value Function Converged, Iteration = \n",
" 3\n",
"\n",
"历时 132.846288 秒。\n",
" 3.728256585425100\n",
"\n",
"Main loop iteration\n",
" 8\n",
"\n",
"Value Function Converged, Iteration = \n",
" 3\n",
"\n",
"历时 132.934235 秒。\n",
" 3.737474350115293\n",
"\n",
"Main loop iteration\n",
" 9\n",
"\n",
"Value Function Converged, Iteration = \n",
" 3\n",
"\n",
"历时 133.043068 秒。\n",
" 3.739684114624203\n",
"\n",
"Main loop iteration\n",
" 10\n",
"\n",
"Value Function Converged, Iteration = \n",
" 2\n",
"\n",
"历时 133.111246 秒。\n",
" 3.740200299833476\n",
"\n",
"Main loop iteration\n",
" 11\n",
"\n",
"Value Function Converged, Iteration = \n",
" 2\n",
"\n",
"历时 133.189599 秒。\n",
" 3.740336725672957\n",
"\n",
"Main loop iteration\n",
" 12\n",
"\n",
"Value Function Converged, Iteration = \n",
" 1\n",
"\n",
"历时 133.240274 秒。\n",
" 3.740342716315000\n",
"\n",
"Main loop iteration\n",
" 13\n",
"\n",
"Value Function Converged, Iteration = \n",
" 1\n",
"\n",
"历时 133.277754 秒。\n",
" 3.740367790362063\n",
"\n"
]
}
],
"source": [
"for iter=1:maxitK\n",
"disp('Main loop iteration')\n",
"disp(iter)\n",
"\n",
" % HAMILTON-JACOBI-BELLMAN EQUATION %\n",
" for n=1:maxit\n",
" V = v;\n",
" % forward difference\n",
" Vaf(1:I-1,:) = (V(2:I,:)-V(1:I-1,:))/da;\n",
" Vaf(I,:) = (w*z + r.*amax).^(-ga); %will never be used, but impose state constraint a<=amax just in case\n",
" % backward difference\n",
" Vab(2:I,:) = (V(2:I,:)-V(1:I-1,:))/da;\n",
" Vab(1,:) = (w*z + r.*amin).^(-ga); %state constraint boundary condition\n",
"\n",
" I_concave = Vab > Vaf; %indicator whether value function is concave (problems arise if this is not the case)\n",
"\n",
" %consumption and savings with forward difference\n",
" cf = Vaf.^(-1/ga);\n",
" sf = w*zz + r.*aa - cf;\n",
" %consumption and savings with backward difference\n",
" cb = Vab.^(-1/ga);\n",
" sb = w*zz + r.*aa - cb;\n",
" %consumption and derivative of value function at steady state\n",
" c0 = w*zz + r.*aa;\n",
" Va0 = c0.^(-ga);\n",
"\n",
" % dV_upwind makes a choice of forward or backward differences based on\n",
" % the sign of the drift\n",
" If = sf > 0; %positive drift --> forward difference\n",
" Ib = sb < 0; %negative drift --> backward difference\n",
" I0 = (1-If-Ib); %at steady state\n",
" %make sure backward difference is used at amax\n",
" % Ib(I,:) = 1; If(I,:) = 0;\n",
" %STATE CONSTRAINT at amin: USE BOUNDARY CONDITION UNLESS sf > 0:\n",
" %already taken care of automatically\n",
"\n",
" Va_Upwind = Vaf.*If + Vab.*Ib + Va0.*I0; %important to include third term\n",
"\n",
" c = Va_Upwind.^(-1/ga);\n",
" u = c.^(1-ga)/(1-ga);\n",
"\n",
" %CONSTRUCT MATRIX A\n",
" X = - min(sb,0)/da;\n",
" Y = - max(sf,0)/da + min(sb,0)/da;\n",
" Z = max(sf,0)/da;\n",
" \n",
" updiag=[0]; %This is needed because of the peculiarity of spdiags.\n",
" for j=1:J\n",
" updiag=[updiag;Z(1:I-1,j);0];\n",
" end\n",
" \n",
" centdiag=reshape(Y,I*J,1);\n",
" \n",
" lowdiag=X(2:I,1);\n",
" for j=2:J\n",
" lowdiag=[lowdiag;0;X(2:I,j)];\n",
" end\n",
" \n",
" AA=spdiags(centdiag,0,I*J,I*J)+spdiags([updiag;0],1,I*J,I*J)+spdiags([lowdiag;0],-1,I*J,I*J);\n",
" \n",
" A = AA + Aswitch;\n",
" \n",
" if max(abs(sum(A,2)))>10^(-9)\n",
" disp('Improper Transition Matrix')\n",
" break\n",
" end\n",
" \n",
" B = (1/Delta + rho)*speye(I*J) - A;\n",
"\n",
" u_stacked = reshape(u,I*J,1);\n",
" V_stacked = reshape(V,I*J,1);\n",
"\n",
" b = u_stacked + V_stacked/Delta;\n",
"\n",
" V_stacked = B\\b; %SOLVE SYSTEM OF EQUATIONS\n",
"\n",
" V = reshape(V_stacked,I,J);\n",
"\n",
" Vchange = V - v;\n",
" v = V;\n",
"\n",
" dist(n) = max(max(abs(Vchange)));\n",
" if dist(n)z_1$.\n",
"\n",
"The process jumps from state 1 to state 2 with intensity $\\lambda_1$ and vice versa with intensity $\\lambda_2$. The two states can be interpreted as employment and unemployment so that $\\lambda_1$ is the job-finding rate and $\\lambda_2$ the job destruction rate."
]
},
{
"cell_type": "markdown",
"id": "bfd5b25a-2a5b-4d70-83fb-a3679b2c70f8",
"metadata": {},
"source": [
"**New change 1: two dimension continuity $\\to$ one continuity plus one discrete**\n",
"\n",
"Compared with O-U process, **Poisson Process is easier**, since productivity/income process now is discrete, so:\n",
"* we don't need to use **Upwind Differential** to analyze how productivity process affect value function.\n",
"* there are only two states, compared to infinite continous O-U process. So, **Grid construction** is different, just use what we do in discrete case.\n",
"* **Transition matrix** is also easier. "
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "7b91ff3b-2fbd-44a7-b29c-19cb4de554e0",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"ans = 'C:\\Users\\ading\\A_TA_Notes'
"
],
"text/plain": [
"ans = 'C:\\Users\\ading\\A_TA_Notes'"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"cd\"C:\\Users\\ading\\A_TA_Notes\"\n",
"pwd"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "131262c3-163b-46b9-b036-cbb53250a100",
"metadata": {},
"outputs": [],
"source": [
"clear all;\n",
"clc; \n",
"close all;\n",
"tic;"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "cef134b2-eac9-41d1-b24c-ac74eec40e39",
"metadata": {},
"outputs": [],
"source": [
"%% Parameters\n",
"\n",
"% households and firms\n",
"ga = 2; % risk aversion\n",
"rho = 0.05; % discount factor\n",
"d = 0.05; % depreciation rate \n",
"al = 1/3; % share for capital\n",
"Aprod = 0.1; % A parameter for production function\n",
"\n",
"% poisson income process\n",
"z1 = 1; % wage state 1\n",
"z2 = 2*z1; % wage state 2\n",
"z = [z1,z2];\n",
"la1 = 1/3; % transfer probability\n",
"la2 = 1/3;\n",
"la = [la1,la2];\n",
"z_ave = (z1*la2 + z2*la1)/(la1 + la2);\n",
"\n",
"%% iteration preparation\n",
"maxit= 100; % max iteration\n",
"crit = 10^(-6); % criterion\n",
"Delta = 1000;\n",
"\n",
"Ir = 40;\n",
"crit_S = 10^(-5);\n",
"\n",
"rmax = 0.049;\n",
"r = 0.04;\n",
"w = 0.05;\n",
"\n",
"r0 = 0.03;\n",
"rmin = 0.01;\n",
"rmax = 0.99*rho;"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "dbe256c6-2e6a-4bde-adc3-50b9cc26dd6e",
"metadata": {},
"outputs": [],
"source": [
"%% set grids\n",
"\n",
"% asset\n",
"I= 1000;\n",
"amin = 0;\n",
"amax = 20;\n",
"a = linspace(amin,amax,I)';\n",
"da = (amax-amin)/(I-1);\n",
"aa = [a,a];\n",
"\n",
"% productivity/income\n",
"z = [z1,z2];\n",
"zz = ones(I,1)*z;\n",
"\n",
"% Preallocation\n",
"% Finite difference approximation of partial derivatives\n",
"% now we only provide this for asset, just one dimension\n",
"% no more approximation for productivity/income\n",
"dVf = zeros(I,2); % derivative value forward\n",
"dVb = zeros(I,2); % derivative value backward\n",
"c = zeros(I,2);"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "370c769c-e137-4186-a2c1-3381d301f794",
"metadata": {},
"outputs": [],
"source": [
"%% Transiton matrix\n",
"Aswitch = [-speye(I)*la(1),speye(I)*la(1);speye(I)*la(2),-speye(I)*la(2)];"
]
},
{
"cell_type": "markdown",
"id": "403561eb-055d-4434-8ce7-0a8eae420a55",
"metadata": {},
"source": [
"**New change 2: Adjusment to HJB equation**"
]
},
{
"cell_type": "markdown",
"id": "aa6d3889-b9e6-45fa-aa5e-045c258bf59c",
"metadata": {},
"source": [
"This change also results from Poisson process."
]
},
{
"cell_type": "markdown",
"id": "6c84199f-2386-4360-ba09-53b554c0dea6",
"metadata": {},
"source": [
"Now the HJB equation is:\n",
"$$\n",
"\\begin{aligned}\n",
"& \\frac{v_{i, j}^{n+1}-v_{i, j}^n}{\\Delta}+\\rho v_{i, j}^{n+1}=u\\left(c_{i, j}^n\\right)+v_{i-1, j}^{n+1} x_{i, j}+v_{i, j}^{n+1} y_{i, j}+v_{i+1, j}^{n+1} z_{i, j}+v_{i,-j}^{n+1} \\lambda_j \\quad \\text { where } \\\\\n",
"& x_{i, j}=-\\frac{\\left(s_{i, j, B}^n\\right)^{-}}{\\Delta a} \\\\\n",
"& y_{i, j}=-\\frac{\\left(s_{i, j, F}^n\\right)^{+}}{\\Delta a}+\\frac{\\left(s_{i, j, B}^n\\right)^{-}}{\\Delta a}-\\lambda_j \\\\\n",
"& z_{i, j}=\\frac{\\left(s_{i, j, F}^n\\right)^{+}}{\\Delta a}\n",
"\\end{aligned}\n",
"$$\n",
"\n",
"Note that importantly $x_{1, j}=z_{I, j}=0, j=1,2$ so $v_{0, j}^{n+1}$ and $v_{I+1, j}^{n+1}$ are never used aboveation (14) is a system of $2 \\times I$ linear equations which can be written in matrix notation as:\n",
"$$\n",
"\\frac{1}{\\Delta}\\left(v^{n+1}-v^n\\right)+\\rho v^{n+1}=u^n+\\mathbf{A}^n v^{n+1}\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "91bf1d8d-b1cf-4cb0-a314-29cbd96c9d83",
"metadata": {},
"source": [
"where $A = x_{i,j} + y_{i,j} + z_{i,j}$"
]
},
{
"cell_type": "markdown",
"id": "c3c0ec3b-d385-4fae-baf7-b59a49d8dc77",
"metadata": {},
"source": [
"Here is the iteration, the code iterate on interest rate."
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "3ab268b8-7864-4d41-8f3a-beeaba8643c9",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Value Function Converged, Iteration = \n",
" 7\n",
"\n",
"历时 2559.207783 秒。\n",
"Excess Demand\n",
"Value Function Converged, Iteration = \n",
" 6\n",
"\n",
"历时 2559.284936 秒。\n",
"Excess Demand\n",
"Value Function Converged, Iteration = \n",
" 5\n",
"\n",
"历时 2559.329921 秒。\n",
"Excess Supply\n",
"Value Function Converged, Iteration = \n",
" 5\n",
"\n",
"历时 2559.365684 秒。\n",
"Excess Supply\n",
"Value Function Converged, Iteration = \n",
" 5\n",
"\n",
"历时 2559.407666 秒。\n",
"Excess Supply\n",
"Value Function Converged, Iteration = \n",
" 5\n",
"\n",
"历时 2559.474438 秒。\n",
"Excess Supply\n",
"Value Function Converged, Iteration = \n",
" 5\n",
"\n",
"历时 2559.512612 秒。\n",
"Excess Demand\n",
"Value Function Converged, Iteration = \n",
" 4\n",
"\n",
"历时 2559.540382 秒。\n",
"Excess Demand\n",
"Value Function Converged, Iteration = \n",
" 4\n",
"\n",
"历时 2559.561689 秒。\n",
"Excess Supply\n",
"Value Function Converged, Iteration = \n",
" 4\n",
"\n",
"历时 2559.587406 秒。\n",
"Excess Demand\n",
"Value Function Converged, Iteration = \n",
" 4\n",
"\n",
"历时 2559.615126 秒。\n",
"Excess Supply\n",
"Value Function Converged, Iteration = \n",
" 4\n",
"\n",
"历时 2559.639559 秒。\n",
"Excess Supply\n",
"Value Function Converged, Iteration = \n",
" 4\n",
"\n",
"历时 2559.671812 秒。\n",
"Excess Supply\n",
"Value Function Converged, Iteration = \n",
" 3\n",
"\n",
"历时 2559.697452 秒。\n",
"Excess Supply\n",
"Value Function Converged, Iteration = \n",
" 3\n",
"\n",
"历时 2559.719571 秒。\n",
"Excess Demand\n",
"Value Function Converged, Iteration = \n",
" 3\n",
"\n",
"历时 2559.741214 秒。\n",
"Equilibrium Found, Interest rate =\n",
" 0.044992080688477\n",
"\n"
]
}
],
"source": [
"%% iteration\n",
"\n",
"for ir=1:Ir;\n",
"\n",
"r_r(ir)=r;\n",
"rmin_r(ir)=rmin;\n",
"rmax_r(ir)=rmax;\n",
"\n",
"KD(ir) = (al*Aprod/(r + d))^(1/(1-al))*z_ave; \n",
"% K(r+d) = alpha * A * K^alpha * Z^(1-alpha)\n",
"w = (1-al)*Aprod*KD(ir).^al*z_ave^(-al);\n",
"% w*Z = (1-alpha) * A * K^alpha * Z^(1-alpha)\n",
"\n",
"if w*z(1) + r*amin < 0\n",
" disp('CAREFUL: borrowing constraint too loose')\n",
"end\n",
"\n",
"v0(:,1) = (w*z(1) + r.*a).^(1-ga)/(1-ga)/rho;\n",
"v0(:,2) = (w*z(2) + r.*a).^(1-ga)/(1-ga)/rho;\n",
"\n",
"if ir>1\n",
"v0 = V_r(:,:,ir-1);\n",
"end\n",
"\n",
"v = v0;\n",
"\n",
"for n=1:maxit\n",
" V = v;\n",
" V_n(:,:,n)=V;\n",
" % forward difference\n",
" dVf(1:I-1,:) = (V(2:I,:)-V(1:I-1,:))/da;\n",
" dVf(I,:) = (w*z + r.*amax).^(-ga); \n",
" %will never be used, but impose state constraint a<=amax just in case\n",
" \n",
" % backward difference\n",
" dVb(2:I,:) = (V(2:I,:)-V(1:I-1,:))/da;\n",
" dVb(1,:) = (w*z + r.*amin).^(-ga); \n",
" %state constraint boundary condition\n",
" \n",
" %consumption and savings with forward difference\n",
" cf = dVf.^(-1/ga);\n",
" ssf = w*zz + r.*aa - cf;\n",
"\n",
" %consumption and savings with backward difference\n",
" cb = dVb.^(-1/ga);\n",
" ssb = w*zz + r.*aa - cb;\n",
"\n",
" %consumption and derivative of value function at steady state\n",
" c0 = w*zz + r.*aa;\n",
" \n",
" % dV_upwind makes a choice of forward or backward differences based on\n",
" % the sign of the drift \n",
" If = ssf > 0; %positive drift --> forward difference\n",
" Ib = ssb < 0; %negative drift --> backward difference\n",
" I0 = (1-If-Ib); %at steady state\n",
" \n",
" c = cf.*If + cb.*Ib + c0.*I0;\n",
" u = c.^(1-ga)/(1-ga);\n",
" \n",
" % Construct Matrix \n",
" X = -min(ssb,0)/da;\n",
" Y = -max(ssf,0)/da + min(ssb,0)/da;\n",
" Z = max(ssf,0)/da;\n",
" \n",
" A1=spdiags(Y(:,1),0,I,I)+spdiags(X(2:I,1),-1,I,I)+spdiags([0;Z(1:I-1,1)],1,I,I);\n",
" A2=spdiags(Y(:,2),0,I,I)+spdiags(X(2:I,2),-1,I,I)+spdiags([0;Z(1:I-1,2)],1,I,I);\n",
" A = [A1,sparse(I,I);sparse(I,I),A2] + Aswitch;\n",
" % S = spdiags( Bin , d , m , n ) creates an m -by- n sparse matrix S \n",
" % by taking the columns of Bin and placing them along the diagonals \n",
" % specified by d \n",
"\n",
" if max(abs(sum(A,2)))>10^(-9)\n",
" disp('Improper Transition Matrix')\n",
" %break\n",
" end\n",
" \n",
" B = (1/Delta + rho)*speye(2*I) - A;\n",
"\n",
" u_stacked = [u(:,1);u(:,2)];\n",
" V_stacked = [V(:,1);V(:,2)];\n",
" \n",
" b = u_stacked + V_stacked/Delta;\n",
" V_stacked = B\\b; % Solve system of equations\n",
" \n",
" V = [V_stacked(1:I), V_stacked(I+1:2*I)];\n",
" \n",
" Vchange = V - v;\n",
" v = V;\n",
"\n",
" dist(n) = max(max(abs(Vchange)));\n",
" if dist(n)crit_S\n",
" disp('Excess Supply')\n",
" rmax = r;\n",
" r = 0.5*(r+rmin);\n",
"elseif S(ir)<-crit_S;\n",
" disp('Excess Demand')\n",
" rmin = r;\n",
" r = 0.5*(r+rmax);\n",
"elseif abs(S(ir))