{ "cells": [ { "cell_type": "markdown", "id": "64cd8102", "metadata": {}, "source": [ "# Estimating an RBC Model with a Stochastic Trend\n", "\n", "Real macro time series have *trends*. US per-capita real GDP, consumption, investment, and the capital stock have all grown at roughly 2% per year for over a century, so the raw data does not fluctuate around a constant mean. The plain RBC model from the previous notebook ([Fitting an RBC to US Data](fit_rbc_to_us_data.ipynb)) has no growth at all — its steady state is a single point, and every endogenous variable is stationary around it. To use that model with real data we had to remove the trend ourselves: log everything, fit a quadratic trend by OLS, and feed the residuals into the Kalman filter. That works, but it has two unattractive features:\n", "\n", "1. **It separates trend and cycle by hand.** Any signal the trend regression absorbs is lost to the structural model — the DSGE never sees it, so it can't tell us anything about it.\n", "2. **It commits us to a trend specification that has nothing to do with economics.** A quadratic-in-time trend is convenient but arbitrary; the data did not generate themselves that way.\n", "\n", "The cleaner alternative is to let the model itself produce the trend. We add a *stochastic trend* to productivity: a permanent random-walk-with-drift component on top of the usual stationary AR(1) shock. The deterministic drift $g_z$ becomes an estimable parameter — interpretable as \"the gross trend growth rate of this economy\" — and trend innovations are real, permanent shocks that show up in the impulse response just like any other shock. After a change-of-variables that divides every trending quantity by the permanent productivity level, the model becomes stationary again and we can solve and estimate it with the same Kalman-filter / log-linearization machinery as the textbook RBC. The estimation is now joint: trend growth and business-cycle parameters are recovered together from the raw (logged) data.\n", "\n", "A trend alone is not enough, though. A bare RBC with a single TFP shock cannot reproduce the joint dynamics of six macro series at once — there is one source of stochastic variation and six things to explain, so the Kalman filter will inevitably hand most of the action to measurement error. To give the model a fighting chance we add three standard real frictions (consumption habit, quadratic investment adjustment cost, variable capital utilisation) and four more structural shocks (preference, labour-disutility, investment, and utilisation-cost shocks). The result is a six-shock, six-observable model: exactly identified before measurement error, with enough flexibility to capture both the trend and the cycle in the French national-accounts data we'll fit it to.\n", "\n", "This notebook walks through the model end-to-end: the equations, the change-of-variables, parser and steady-state checks, then estimation with PyMC." ] }, { "cell_type": "markdown", "id": "f9d0a217", "metadata": {}, "source": [ "
gEconpy. Install it into the same environment before running the cells below:\n",
"\n",
"pip install eurostat\n", "\n", "
eurostat pulls quarterly French national-accounts and labour-force series directly from Eurostat’s SDMX API.\n",
"\\[u_{t} = - \\frac{L_{t}^{\\sigma_{L} + 1} \\Theta_{t}}{\\sigma_{L} + 1} + \\log{\\left(- \\frac{\\phi_{H} C_{t-1}}{G_{t}} + C_{t} \\right)}\\]
\n", "\\[\\Psi_{z t} = shock_{\\psi t} \\left(\\frac{\\psi_{z} \\left(z_{t} - 1\\right)^{2}}{2} + \\psi_{z 1} \\left(z_{t} - 1\\right)\\right)\\]
\n", "\\[adj_{cost t} = - \\frac{\\gamma_{I} \\left(g_{z} - \\frac{G_{t} I_{t} shock_{I t}}{I_{t-1}}\\right)^{2}}{2} + 1\\]
\n", "\\[\\operatorname{Set}\\left(\\left[ C_{t}, \\ L_{t}, \\ I_{t}, \\ K_{t}, \\ z_{t}\\right]\\right)\\]
\n", "\\[U_{t} = U_{t+1} \\beta_{t} + u_{t}\\]
\n", "\\[C_{t} + I_{t} + \\frac{K_{t-1} \\Psi_{z t}}{G_{t}} = L_{t} W_{t} + \\frac{K_{t-1} r_{t} z_{t}}{G_{t}}\\]
\n", "\\[K_{t} = I_{t} adj_{cost t} - \\frac{K_{t-1} \\left(\\delta - 1\\right)}{G_{t}}\\]
\n", "\\[\\log{\\left(\\beta_{t} \\right)} = \\rho_{\\beta} \\log{\\left(\\beta_{t-1} \\right)} + \\epsilon_{\\beta t} - \\left(\\rho_{\\beta} - 1\\right) \\log{\\left(\\beta \\right)}\\]
\n", "\\[\\log{\\left(\\Theta_{t} \\right)} = \\rho_{\\Theta} \\log{\\left(\\Theta_{t-1} \\right)} + \\epsilon_{\\Theta t} - \\left(\\rho_{\\Theta} - 1\\right) \\log{\\left(\\Theta \\right)}\\]
\n", "\\[\\log{\\left(shock_{I t} \\right)} = \\rho_{I} \\log{\\left(shock_{I t-1} \\right)} + \\epsilon_{I t}\\]
\n", "\\[\\log{\\left(shock_{\\psi t} \\right)} = \\rho_{\\psi} \\log{\\left(shock_{\\psi t-1} \\right)} + \\epsilon_{\\psi t}\\]
\n", "\\[R_{f t} + 1 = \\frac{G_{t+1} \\lambda_{t}}{\\beta_{t} \\lambda_{t+1}}\\]
\n", "\\[\\operatorname{Set}\\left(\\left[ \\epsilon_{\\beta t}, \\ \\epsilon_{\\Theta t}, \\ \\epsilon_{I t}, \\ \\epsilon_{\\psi t}\\right]\\right)\\]
\n", "\\[\\beta = 0.99\\]
\n", "\\[\\delta = 0.025\\]
\n", "\\[\\sigma_{L} = 2.0\\]
\n", "\\[\\phi_{H} = 0.6\\]
\n", "\\[\\gamma_{I} = 6.0\\]
\n", "\\[\\psi_{z} = 0.2\\]
\n", "\\[\\psi_{z 1} = \\delta - 1 + \\frac{g_{z}}{\\beta}\\]
\n", "\\[\\Theta = 1\\]
\n", "\\[\\rho_{\\beta} = 0.9\\]
\n", "\\[\\rho_{\\Theta} = 0.9\\]
\n", "\\[\\rho_{I} = 0.9\\]
\n", "\\[\\rho_{\\psi} = 0.9\\]
\n", "\\[\\operatorname{Set}\\left(\\left[ K_{d t}, \\ L_{t}\\right]\\right)\\]
\n", "\\[TC_{t} = - K_{d t} r_{t} - L_{t} W_{t}\\]
\n", "\\[Y_{t} = A_{t} K_{d t}^{\\alpha} L_{t}^{1 - \\alpha}\\]
\n", "\\[mc_{t} = 1\\]
\n", "\\[K_{d t} = \\frac{K_{t-1} z_{t}}{G_{t}}\\]
\n", "\\[\\alpha = 0.3\\]
\n", "\\[\\log{\\left(A_{t} \\right)} = \\rho_{A} \\log{\\left(A_{t-1} \\right)} + \\epsilon_{A t}\\]
\n", "\\[\\log{\\left(G_{t} \\right)} = \\epsilon_{z t} + \\log{\\left(g_{z} \\right)}\\]
\n", "\\[\\operatorname{Set}\\left(\\left[ \\epsilon_{A t}, \\ \\epsilon_{z t}\\right]\\right)\\]
\n", "\\[\\rho_{A} = 0.95\\]
\n", "\\[g_{z} = 1.005\\]
\n", "| \n", " | value | \n", "
|---|---|
| A_ss | \n", "1.0000 | \n", "
| C_ss | \n", "1.7838 | \n", "
| G_ss | \n", "1.0050 | \n", "
| I_ss | \n", "0.5153 | \n", "
| K_ss | \n", "17.2641 | \n", "
| L_ss | \n", "0.9710 | \n", "
| R_f_ss | \n", "0.0152 | \n", "
| Theta_ss | \n", "1.0000 | \n", "
| W_ss | \n", "1.6574 | \n", "
| Y_ss | \n", "2.2991 | \n", "
| beta_ss | \n", "0.9900 | \n", "
| lambda_ss | \n", "0.5689 | \n", "
| q_ss | \n", "0.5689 | \n", "
| r_ss | \n", "0.0402 | \n", "
| shock_I_ss | \n", "1.0000 | \n", "
| shock_psi_ss | \n", "1.0000 | \n", "
| z_ss | \n", "1.0000 | \n", "
| \n", " | mean | \n", "sd | \n", "eti89_lb | \n", "eti89_ub | \n", "ess_bulk | \n", "ess_tail | \n", "mcse_mean | \n", "mcse_sd | \n", "
|---|---|---|---|---|---|---|---|---|
| C_ss | \n", "0.6 | \n", "0.2 | \n", "0.29 | \n", "0.94 | \n", "522 | \n", "483 | \n", "0.0086 | \n", "0.0058 | \n", "
| G_ss | \n", "1.005 | \n", "0.002 | \n", "1 | \n", "1 | \n", "621 | \n", "466 | \n", "7.9e-05 | \n", "6e-05 | \n", "
| I_ss | \n", "0.127 | \n", "0.063 | \n", "0.048 | \n", "0.23 | \n", "537 | \n", "462 | \n", "0.0027 | \n", "0.0025 | \n", "
| K_ss | \n", "4.2 | \n", "2.5 | \n", "1.6 | \n", "8.6 | \n", "497 | \n", "433 | \n", "0.12 | \n", "0.16 | \n", "
| L_ss | \n", "0.377 | \n", "0.112 | \n", "0.2 | \n", "0.55 | \n", "542 | \n", "485 | \n", "0.0048 | \n", "0.0028 | \n", "
| R_f_ss | \n", "0.025 | \n", "0.011 | \n", "0.011 | \n", "0.046 | \n", "482 | \n", "367 | \n", "0.0005 | \n", "0.0005 | \n", "
| Theta_ss | \n", "157 | \n", "19 | \n", "130 | \n", "190 | \n", "511 | \n", "504 | \n", "0.86 | \n", "0.58 | \n", "
| W_ss | \n", "1.39 | \n", "0.21 | \n", "1.1 | \n", "1.7 | \n", "499 | \n", "430 | \n", "0.0093 | \n", "0.0092 | \n", "
| Y_ss | \n", "0.72 | \n", "0.25 | \n", "0.34 | \n", "1.1 | \n", "544 | \n", "501 | \n", "0.011 | \n", "0.0078 | \n", "
| beta_ss | \n", "0.9804 | \n", "0.0105 | \n", "0.96 | \n", "0.99 | \n", "492 | \n", "354 | \n", "0.00047 | \n", "0.00044 | \n", "
| lambda_ss | \n", "2.05 | \n", "0.9 | \n", "1.1 | \n", "3.7 | \n", "531 | \n", "438 | \n", "0.04 | \n", "0.053 | \n", "
| q_ss | \n", "2.05 | \n", "0.9 | \n", "1.1 | \n", "3.7 | \n", "531 | \n", "438 | \n", "0.04 | \n", "0.053 | \n", "
| r_ss | \n", "0.052 | \n", "0.014 | \n", "0.033 | \n", "0.078 | \n", "426 | \n", "382 | \n", "0.00067 | \n", "0.00061 | \n", "
| \n", " | Y | \n", "C | \n", "I | \n", "L_hours | \n", "D1 | \n", "P | \n", "N | \n", "i_nom | \n", "
|---|---|---|---|---|---|---|---|---|
| Time | \n", "\n", " | \n", " | \n", " | \n", " | \n", " | \n", " | \n", " | \n", " |
| 2025-01-01 | \n", "615341.9 | \n", "329412.8 | \n", "133012.2 | \n", "11555800.0 | \n", "381012.9 | \n", "120.120 | \n", "41354.5 | \n", "3.30500 | \n", "
| 2025-04-01 | \n", "617501.9 | \n", "329332.3 | \n", "133526.0 | \n", "11543600.0 | \n", "382597.8 | \n", "120.174 | \n", "41371.6 | \n", "3.25341 | \n", "
| 2025-07-01 | \n", "620923.5 | \n", "329660.8 | \n", "134567.3 | \n", "11550900.0 | \n", "384013.7 | \n", "120.565 | \n", "41392.8 | \n", "3.43000 | \n", "
| 2025-10-01 | \n", "622228.7 | \n", "330832.7 | \n", "134988.3 | \n", "11542600.0 | \n", "386328.9 | \n", "121.106 | \n", "41414.4 | \n", "3.44000 | \n", "
| 2026-01-01 | \n", "622197.9 | \n", "330450.4 | \n", "134487.8 | \n", "NaN | \n", "NaN | \n", "NaN | \n", "NaN | \n", "3.52744 | \n", "
| \n", " | Y | \n", "C | \n", "I | \n", "L | \n", "w | \n", "r | \n", "
|---|---|---|---|---|---|---|
| Time | \n", "\n", " | \n", " | \n", " | \n", " | \n", " | \n", " |
| 2025-01-01 | \n", "615341.9 | \n", "329412.8 | \n", "133012.2 | \n", "0.194050 | \n", "0.027449 | \n", "0.004561 | \n", "
| 2025-04-01 | \n", "617501.9 | \n", "329332.3 | \n", "133526.0 | \n", "0.191636 | \n", "0.027580 | \n", "0.007583 | \n", "
| 2025-07-01 | \n", "620923.5 | \n", "329660.8 | \n", "134567.3 | \n", "0.189576 | \n", "0.027575 | \n", "0.005196 | \n", "
| 2025-10-01 | \n", "622228.7 | \n", "330832.7 | \n", "134988.3 | \n", "0.189341 | \n", "0.027637 | \n", "0.003986 | \n", "
| 2026-01-01 | \n", "622197.9 | \n", "330450.4 | \n", "134487.8 | \n", "NaN | \n", "NaN | \n", "NaN | \n", "
| \n", " | dlog_Y | \n", "dlog_C | \n", "dlog_I | \n", "dlog_W | \n", "log_L | \n", "R_f | \n", "
|---|---|---|---|---|---|---|
| count | \n", "184.000000 | \n", "184.000000 | \n", "184.000000 | \n", "183.000000 | \n", "97.000000 | \n", "183.000000 | \n", "
| mean | \n", "0.004174 | \n", "0.004232 | \n", "0.003785 | \n", "0.003108 | \n", "-1.721223 | \n", "0.007613 | \n", "
| std | \n", "0.015697 | \n", "0.016851 | \n", "0.019370 | \n", "0.008499 | \n", "0.039445 | \n", "0.010591 | \n", "
| min | \n", "-0.129838 | \n", "-0.108797 | \n", "-0.134123 | \n", "-0.049190 | \n", "-1.934368 | \n", "-0.033394 | \n", "
| 25% | \n", "0.001593 | \n", "0.000756 | \n", "-0.002873 | \n", "0.000461 | \n", "-1.742982 | \n", "0.001986 | \n", "
| 50% | \n", "0.004425 | \n", "0.004226 | \n", "0.004512 | \n", "0.003426 | \n", "-1.728491 | \n", "0.006082 | \n", "
| 75% | \n", "0.007129 | \n", "0.007909 | \n", "0.011448 | \n", "0.005755 | \n", "-1.709883 | \n", "0.011844 | \n", "
| max | \n", "0.142350 | \n", "0.158507 | \n", "0.152297 | \n", "0.064056 | \n", "-1.639637 | \n", "0.066406 | \n", "
| \n", " | dlog_Y | \n", "dlog_C | \n", "dlog_I | \n", "dlog_W | \n", "log_L | \n", "R_f | \n", "
|---|---|---|---|---|---|---|
| Time | \n", "\n", " | \n", " | \n", " | \n", " | \n", " | \n", " |
| 2019-10-01 | \n", "-0.005059 | \n", "-0.00154 | \n", "-0.000271 | \n", "-0.002846 | \n", "-1.706902 | \n", "-0.003547 | \n", "
| 2020-01-01 | \n", "NaN | \n", "NaN | \n", "NaN | \n", "NaN | \n", "NaN | \n", "NaN | \n", "
| 2020-04-01 | \n", "NaN | \n", "NaN | \n", "NaN | \n", "NaN | \n", "NaN | \n", "NaN | \n", "
| 2020-07-01 | \n", "NaN | \n", "NaN | \n", "NaN | \n", "NaN | \n", "NaN | \n", "NaN | \n", "
| 2020-10-01 | \n", "NaN | \n", "NaN | \n", "NaN | \n", "NaN | \n", "NaN | \n", "NaN | \n", "
| 2021-01-01 | \n", "NaN | \n", "NaN | \n", "NaN | \n", "NaN | \n", "NaN | \n", "NaN | \n", "
| 2021-04-01 | \n", "NaN | \n", "NaN | \n", "NaN | \n", "NaN | \n", "NaN | \n", "NaN | \n", "
| 2021-07-01 | \n", "NaN | \n", "NaN | \n", "NaN | \n", "NaN | \n", "NaN | \n", "NaN | \n", "
Model Requirements \n", " \n", " Variable Shape Constraints Dimensions \n", " ──────────────────────────────────────────────────────── \n", " Theta () None \n", " alpha () None \n", " beta () None \n", " delta () None \n", " g_z () None \n", " gamma_I () None \n", " phi_H () None \n", " psi_z () None \n", " rho_A () None \n", " rho_I () None \n", " rho_Theta () None \n", " rho_beta () None \n", " rho_psi () None \n", " sigma_L () None \n", " sigma_epsilon_A () Positive None \n", " sigma_epsilon_I () Positive None \n", " sigma_epsilon_Theta () Positive None \n", " sigma_epsilon_beta () Positive None \n", " sigma_epsilon_psi () Positive None \n", " sigma_epsilon_z () Positive None \n", " error_sigma_dlog_Y () Positive None \n", " error_sigma_dlog_C () Positive None \n", " error_sigma_dlog_I () Positive None \n", " error_sigma_dlog_W () Positive None \n", " error_sigma_log_L () Positive None \n", " error_sigma_R_f () Positive None \n", " \n", " These parameters should be assigned priors inside a PyMC \n", " model block before calling the build_statespace_graph \n", " method. \n", "\n" ], "text/plain": [ "\u001b[3m Model Requirements \u001b[0m\n", " \n", " \u001b[1m \u001b[0m\u001b[1mVariable \u001b[0m\u001b[1m \u001b[0m \u001b[1m \u001b[0m\u001b[1mShape\u001b[0m\u001b[1m \u001b[0m \u001b[1m \u001b[0m\u001b[1mConstraints\u001b[0m\u001b[1m \u001b[0m \u001b[1m \u001b[0m\u001b[1mDimensions\u001b[0m\u001b[1m \u001b[0m \n", " ──────────────────────────────────────────────────────── \n", " Theta \u001b[1m(\u001b[0m\u001b[1m)\u001b[0m \u001b[3;35mNone\u001b[0m \n", " alpha \u001b[1m(\u001b[0m\u001b[1m)\u001b[0m \u001b[3;35mNone\u001b[0m \n", " beta \u001b[1m(\u001b[0m\u001b[1m)\u001b[0m \u001b[3;35mNone\u001b[0m \n", " delta \u001b[1m(\u001b[0m\u001b[1m)\u001b[0m \u001b[3;35mNone\u001b[0m \n", " g_z \u001b[1m(\u001b[0m\u001b[1m)\u001b[0m \u001b[3;35mNone\u001b[0m \n", " gamma_I \u001b[1m(\u001b[0m\u001b[1m)\u001b[0m \u001b[3;35mNone\u001b[0m \n", " phi_H \u001b[1m(\u001b[0m\u001b[1m)\u001b[0m \u001b[3;35mNone\u001b[0m \n", " psi_z \u001b[1m(\u001b[0m\u001b[1m)\u001b[0m \u001b[3;35mNone\u001b[0m \n", " rho_A \u001b[1m(\u001b[0m\u001b[1m)\u001b[0m \u001b[3;35mNone\u001b[0m \n", " rho_I \u001b[1m(\u001b[0m\u001b[1m)\u001b[0m \u001b[3;35mNone\u001b[0m \n", " rho_Theta \u001b[1m(\u001b[0m\u001b[1m)\u001b[0m \u001b[3;35mNone\u001b[0m \n", " rho_beta \u001b[1m(\u001b[0m\u001b[1m)\u001b[0m \u001b[3;35mNone\u001b[0m \n", " rho_psi \u001b[1m(\u001b[0m\u001b[1m)\u001b[0m \u001b[3;35mNone\u001b[0m \n", " sigma_L \u001b[1m(\u001b[0m\u001b[1m)\u001b[0m \u001b[3;35mNone\u001b[0m \n", " sigma_epsilon_A \u001b[1m(\u001b[0m\u001b[1m)\u001b[0m Positive \u001b[3;35mNone\u001b[0m \n", " sigma_epsilon_I \u001b[1m(\u001b[0m\u001b[1m)\u001b[0m Positive \u001b[3;35mNone\u001b[0m \n", " sigma_epsilon_Theta \u001b[1m(\u001b[0m\u001b[1m)\u001b[0m Positive \u001b[3;35mNone\u001b[0m \n", " sigma_epsilon_beta \u001b[1m(\u001b[0m\u001b[1m)\u001b[0m Positive \u001b[3;35mNone\u001b[0m \n", " sigma_epsilon_psi \u001b[1m(\u001b[0m\u001b[1m)\u001b[0m Positive \u001b[3;35mNone\u001b[0m \n", " sigma_epsilon_z \u001b[1m(\u001b[0m\u001b[1m)\u001b[0m Positive \u001b[3;35mNone\u001b[0m \n", " error_sigma_dlog_Y \u001b[1m(\u001b[0m\u001b[1m)\u001b[0m Positive \u001b[3;35mNone\u001b[0m \n", " error_sigma_dlog_C \u001b[1m(\u001b[0m\u001b[1m)\u001b[0m Positive \u001b[3;35mNone\u001b[0m \n", " error_sigma_dlog_I \u001b[1m(\u001b[0m\u001b[1m)\u001b[0m Positive \u001b[3;35mNone\u001b[0m \n", " error_sigma_dlog_W \u001b[1m(\u001b[0m\u001b[1m)\u001b[0m Positive \u001b[3;35mNone\u001b[0m \n", " error_sigma_log_L \u001b[1m(\u001b[0m\u001b[1m)\u001b[0m Positive \u001b[3;35mNone\u001b[0m \n", " error_sigma_R_f \u001b[1m(\u001b[0m\u001b[1m)\u001b[0m Positive \u001b[3;35mNone\u001b[0m \n", " \n", "\u001b[2;3m These parameters should be assigned priors inside a PyMC \u001b[0m\n", "\u001b[2;3m model block before calling the build_statespace_graph \u001b[0m\n", "\u001b[2;3m method. \u001b[0m\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "observation_equations = {\n", " \"dlog_Y\": \"log(Y[] / Y[-1] * G[])\",\n", " \"dlog_C\": \"log(C[] / C[-1] * G[])\",\n", " \"dlog_I\": \"log(I[] / I[-1] * G[])\",\n", " \"dlog_W\": \"log(W[] / W[-1] * G[])\",\n", " \"log_L\": \"log(L[])\",\n", "}\n", "\n", "# R_f has no observation equation: it is observed directly off the latent\n", "# state, with the steady-state value used as the observation intercept.\n", "ss_mod.configure(\n", " observed_states=list(df_obs.columns),\n", " observation_equations=observation_equations,\n", " ss_obs_intercept=[\"R_f\"],\n", " measurement_error=list(df_obs.columns),\n", " mode=\"JAX\",\n", " solver=\"scan_cycle_reduction\",\n", " max_iter=13,\n", " use_adjoint_gradients=True,\n", " use_direct_lyapunov=False,\n", ")" ] }, { "cell_type": "markdown", "id": "1dbaea72", "metadata": {}, "source": [ "### Building the PyMC graph\n", "\n", "Every parameter listed by `configure()` needs to live as a PyMC random variable before we can call `build_statespace_graph`. The 14 structural parameters already have priors from the GCN, so `ss_mod.to_pymc()` registers them in one call. The remaining twelve standard-deviations (six structural shocks, six measurement errors) are not in the GCN by design — they're estimation-time choices — so we set them by hand.\n", "\n", "The data is in growth-rate / log-level scale, so both kinds of $\\sigma$ live somewhere around the 0.001–0.05 band; a `maxent(Gamma, lower=0.001, upper=0.05)` prior pins 99% of the prior mass inside that range." ] }, { "cell_type": "code", "execution_count": 19, "id": "8b9eec12", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/jessegrabowski/mambaforge/envs/geconpy-dev/lib/python3.12/site-packages/pymc_extras/statespace/utils/data_tools.py:179: ImputationWarning: Provided data contains missing values and will be automatically imputed as hidden states during Kalman filtering.\n", " warnings.warn(impute_message, ImputationWarning)\n" ] } ], "source": [ "with pm.Model(coords=ss_mod.coords) as rbc_growth_pm:\n", " ss_mod.to_pymc()\n", " prior_dict = ss_mod.param_priors.copy()\n", "\n", " # Theta sets the steady-state level of labour and has no GCN prior. It is\n", " # only weakly identified: L_ss is a joint function of Theta, sigma_L and\n", " # alpha, so a loose prior leaves a ridge the sampler cannot traverse. Pin\n", " # it with a tight Normal centred on the value implied by the observed\n", " # labour mean -- Theta ~ 157 gives L_ss = exp(mean(log_L)) ~ 0.18.\n", " d = pz.Normal(mu=157.0, sigma=10.0)\n", " prior_dict[\"Theta\"] = d\n", " d.to_pymc(\"Theta\")\n", "\n", " for shock in ss_mod.shock_names:\n", " d = pz.maxent(pz.Gamma(), lower=0.001, upper=0.05, plot=False)\n", " prior_dict[f\"sigma_{shock}\"] = d\n", " d.to_pymc(f\"sigma_{shock}\")\n", "\n", " for obs in ss_mod.error_states:\n", " d = pz.maxent(pz.Gamma(), lower=0.001, upper=0.05, plot=False)\n", " prior_dict[f\"error_sigma_{obs}\"] = d\n", " d.to_pymc(f\"error_sigma_{obs}\")\n", "\n", " ss_mod.build_statespace_graph(\n", " data=df_obs_masked[list(ss_mod.observed_states)],\n", " )" ] }, { "cell_type": "code", "execution_count": 20, "id": "c15e399b-325c-47da-89ba-0cd12e289388", "metadata": {}, "outputs": [], "source": [ "import nutpie as ntp\n", "\n", "initial_points = {k: v + 1e-7 for k, v in ss_mod.param_dict.items() if k in [x.name for x in rbc_growth_pm.basic_RVs]}\n", "initial_points[\"Theta\"] = 100.0\n", "ntp_mod = ntp.compile_pymc_model(\n", " rbc_growth_pm,\n", " backend=\"jax\",\n", " gradient_backend=\"jax\",\n", " freeze_model=True,\n", " default_initialization_strategy=\"prior\",\n", " jitter_rvs=None,\n", " initial_points=initial_points,\n", ")" ] }, { "cell_type": "code", "execution_count": 21, "id": "b33c40df-22e8-4d2e-9f3b-cc4ffa3db893", "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "
Sampler Progress
\n", "Total Chains: 4
\n", "Active Chains: 0
\n", "\n", " Finished Chains:\n", " 4\n", "
\n", "Sampling for 7 minutes
\n", "\n", " Estimated Time to Completion:\n", " now\n", "
\n", "\n", " \n", "| Progress | \n", "Draws | \n", "Divergences | \n", "Step Size | \n", "Gradients/Draw | \n", "
|---|---|---|---|---|
| \n", " \n", " | \n", "2500 | \n", "0 | \n", "0.65 | \n", "7 | \n", "
| \n", " \n", " | \n", "2500 | \n", "0 | \n", "0.69 | \n", "7 | \n", "
| \n", " \n", " | \n", "2500 | \n", "0 | \n", "0.69 | \n", "7 | \n", "
| \n", " \n", " | \n", "2500 | \n", "0 | \n", "0.69 | \n", "7 | \n", "
| \n", " | mean | \n", "sd | \n", "eti89_lb | \n", "eti89_ub | \n", "ess_bulk | \n", "ess_tail | \n", "r_hat | \n", "mcse_mean | \n", "mcse_sd | \n", "
|---|---|---|---|---|---|---|---|---|---|
| Theta | \n", "149.908 | \n", "10.207 | \n", "133.412 | \n", "166.172 | \n", "5894.716 | \n", "1588.999 | \n", "1.004 | \n", "0.133 | \n", "0.089 | \n", "
| alpha | \n", "0.215 | \n", "0.021 | \n", "0.183 | \n", "0.249 | \n", "4678.697 | \n", "1391.190 | \n", "1.004 | \n", "0.000 | \n", "0.000 | \n", "
| beta | \n", "0.994 | \n", "0.002 | \n", "0.991 | \n", "0.996 | \n", "4113.808 | \n", "1589.237 | \n", "1.001 | \n", "0.000 | \n", "0.000 | \n", "
| delta | \n", "0.037 | \n", "0.010 | \n", "0.023 | \n", "0.054 | \n", "4665.378 | \n", "1653.002 | \n", "1.005 | \n", "0.000 | \n", "0.000 | \n", "
| g_z | \n", "1.004 | \n", "0.000 | \n", "1.003 | \n", "1.004 | \n", "6545.150 | \n", "1672.232 | \n", "1.005 | \n", "0.000 | \n", "0.000 | \n", "
| gamma_I | \n", "7.887 | \n", "1.031 | \n", "6.312 | \n", "9.563 | \n", "5049.032 | \n", "1624.744 | \n", "1.004 | \n", "0.015 | \n", "0.011 | \n", "
| phi_H | \n", "0.821 | \n", "0.070 | \n", "0.701 | \n", "0.918 | \n", "4542.009 | \n", "1587.434 | \n", "1.000 | \n", "0.001 | \n", "0.001 | \n", "
| psi_z | \n", "0.255 | \n", "0.088 | \n", "0.126 | \n", "0.401 | \n", "4113.406 | \n", "1419.538 | \n", "1.000 | \n", "0.001 | \n", "0.001 | \n", "
| rho_A | \n", "0.931 | \n", "0.035 | \n", "0.870 | \n", "0.976 | \n", "4303.360 | \n", "1511.364 | \n", "1.005 | \n", "0.001 | \n", "0.000 | \n", "
| rho_I | \n", "0.775 | \n", "0.077 | \n", "0.643 | \n", "0.886 | \n", "4458.025 | \n", "1771.862 | \n", "1.004 | \n", "0.001 | \n", "0.001 | \n", "
| rho_Theta | \n", "0.963 | \n", "0.011 | \n", "0.946 | \n", "0.979 | \n", "5820.256 | \n", "1373.250 | \n", "1.008 | \n", "0.000 | \n", "0.000 | \n", "
| rho_beta | \n", "0.595 | \n", "0.087 | \n", "0.461 | \n", "0.732 | \n", "3547.959 | \n", "1528.331 | \n", "1.007 | \n", "0.001 | \n", "0.001 | \n", "
| rho_psi | \n", "0.759 | \n", "0.092 | \n", "0.595 | \n", "0.892 | \n", "5505.926 | \n", "1358.054 | \n", "1.005 | \n", "0.001 | \n", "0.001 | \n", "
| sigma_L | \n", "1.866 | \n", "0.055 | \n", "1.774 | \n", "1.946 | \n", "4792.025 | \n", "1502.964 | \n", "1.001 | \n", "0.001 | \n", "0.001 | \n", "
| sigma_epsilon_A | \n", "0.002 | \n", "0.001 | \n", "0.000 | \n", "0.003 | \n", "4600.018 | \n", "1502.088 | \n", "1.001 | \n", "0.000 | \n", "0.000 | \n", "
| sigma_epsilon_I | \n", "0.011 | \n", "0.002 | \n", "0.007 | \n", "0.015 | \n", "4231.803 | \n", "1677.051 | \n", "1.002 | \n", "0.000 | \n", "0.000 | \n", "
| sigma_epsilon_Theta | \n", "0.024 | \n", "0.003 | \n", "0.020 | \n", "0.029 | \n", "5033.407 | \n", "1497.048 | \n", "1.004 | \n", "0.000 | \n", "0.000 | \n", "
| sigma_epsilon_beta | \n", "0.005 | \n", "0.001 | \n", "0.003 | \n", "0.007 | \n", "3588.853 | \n", "1746.742 | \n", "1.000 | \n", "0.000 | \n", "0.000 | \n", "
| sigma_epsilon_psi | \n", "0.015 | \n", "0.010 | \n", "0.003 | \n", "0.034 | \n", "5285.584 | \n", "1607.809 | \n", "1.001 | \n", "0.000 | \n", "0.000 | \n", "
| sigma_epsilon_z | \n", "0.004 | \n", "0.001 | \n", "0.003 | \n", "0.005 | \n", "3856.080 | \n", "1553.218 | \n", "1.002 | \n", "0.000 | \n", "0.000 | \n", "
| error_sigma_dlog_Y | \n", "0.003 | \n", "0.000 | \n", "0.002 | \n", "0.003 | \n", "4398.187 | \n", "1393.186 | \n", "1.004 | \n", "0.000 | \n", "0.000 | \n", "
| error_sigma_dlog_C | \n", "0.005 | \n", "0.000 | \n", "0.004 | \n", "0.005 | \n", "5222.385 | \n", "1468.102 | \n", "1.003 | \n", "0.000 | \n", "0.000 | \n", "
| error_sigma_dlog_I | \n", "0.004 | \n", "0.001 | \n", "0.002 | \n", "0.005 | \n", "4749.154 | \n", "1370.459 | \n", "1.004 | \n", "0.000 | \n", "0.000 | \n", "
| error_sigma_dlog_W | \n", "0.004 | \n", "0.000 | \n", "0.004 | \n", "0.005 | \n", "5861.110 | \n", "1487.360 | \n", "1.000 | \n", "0.000 | \n", "0.000 | \n", "
| error_sigma_log_L | \n", "0.009 | \n", "0.001 | \n", "0.008 | \n", "0.011 | \n", "6602.060 | \n", "1772.543 | \n", "1.003 | \n", "0.000 | \n", "0.000 | \n", "
| error_sigma_R_f | \n", "0.007 | \n", "0.001 | \n", "0.006 | \n", "0.009 | \n", "4546.410 | \n", "1354.833 | \n", "1.002 | \n", "0.000 | \n", "0.000 | \n", "