{ "cells": [ { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "┌ Info: Precompiling Turing [fce5fe82-541a-59a6-adf8-730c64b5f9a0]\n", "└ @ Base loading.jl:1342\n" ] }, { "data": { "text/plain": [ "21-element Vector{Vector{Float64}}:\n", " [7.824904812740593e-5, 0.0]\n", " [0.009388280821124171, -7.824904812740593e-5]\n", " [0.004141003230806073, 0.08564547473108383]\n", " [0.003113377412500409, 0.043761814038505646]\n", " [0.007700917546593067, 0.03114748946994769]\n", " [0.0028219400485847866, 0.07243491327720691]\n", " [0.014351652682013355, 0.030627394848069485]\n", " [0.014965386945969147, 0.1336674064514162]\n", " [0.014501839251820175, 0.1470204417729682]\n", " [0.012121007622016115, 0.14375484164268565]\n", " [0.006158082144511833, 0.12158897809077801]\n", " [0.0022856243818403987, 0.06503713781067733]\n", " [-0.005567886073177921, 0.02513421406357805]\n", " [-0.007159495107005886, -0.05012136636762045]\n", " [-0.0009141646711710215, -0.07044723363875496]\n", " [0.006966287246748618, -0.014495428593979284]\n", " [0.013428142860868263, 0.06225784779844589]\n", " [0.013777970823039514, 0.1275144960461263]\n", " [0.009452327892967818, 0.13562318681100563]\n", " [0.009420130597109135, 0.09641143723837649]\n", " [0.0071400525559970375, 0.09316115727178581]" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using DifferentiableStateSpaceModels, LinearAlgebra, Turing, Zygote\n", "using DifferentiableStateSpaceModels.Examples\n", "using Turing: @addlogprob!\n", "Turing.setadbackend(:zygote)\n", "\n", "# Create models from modules and then solve\n", "model_rbc = @include_example_module(Examples.rbc_observables)\n", "\n", "# Generate artificial data for estimation\n", "p_f = (ρ=0.2, δ=0.02, σ=0.01, Ω_1=0.01) # Fixed parameters\n", "p_d = (α=0.5, β=0.95) # Pseudo-true values\n", "sol = generate_perturbation(model_rbc, p_d, p_f, Val(1))\n", "sol_second = generate_perturbation(model_rbc, p_d, p_f, Val(2))\n", "\n", "T = 20\n", "ϵ = [randn(model_rbc.n_ϵ) for _ in 1:T]\n", "x0 = zeros(model_rbc.n_x)\n", "fake_z = solve(sol, x0, (0, T), DifferentiableStateSpaceModels.LTI(); noise = ϵ).z\n", "fake_z_second = solve(sol_second, x0, (0, T), DifferentiableStateSpaceModels.QTI(); noise = ϵ).z" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "┌ Info: Found initial step size\n", "│ ϵ = 0.2\n", "└ @ Turing.Inference C:\\Users\\wupei\\.julia\\packages\\Turing\\nfMhU\\src\\inference\\hmc.jl:188\n", "┌ Warning: The current proposal will be rejected due to numerical error(s).\n", "│ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)\n", "└ @ AdvancedHMC C:\\Users\\wupei\\.julia\\packages\\AdvancedHMC\\HQHnm\\src\\hamiltonian.jl:47\n", "┌ Warning: The current proposal will be rejected due to numerical error(s).\n", "│ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)\n", "└ @ AdvancedHMC C:\\Users\\wupei\\.julia\\packages\\AdvancedHMC\\HQHnm\\src\\hamiltonian.jl:47\n", "┌ Warning: The current proposal will be rejected due to numerical error(s).\n", "│ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)\n", "└ @ AdvancedHMC C:\\Users\\wupei\\.julia\\packages\\AdvancedHMC\\HQHnm\\src\\hamiltonian.jl:47\n", "\u001b[32mSampling: 100%|█████████████████████████████████████████| Time: 0:00:42\u001b[39m\n" ] }, { "data": { "text/plain": [ "Chains MCMC chain (1000×14×1 Array{Float64, 3}):\n", "\n", "Iterations = 101:1:1100\n", "Number of chains = 1\n", "Samples per chain = 1000\n", "Wall duration = 121.72 seconds\n", "Compute duration = 121.72 seconds\n", "parameters = α, β\n", "internals = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size\n", "\n", "Summary Statistics\n", " \u001b[1m parameters \u001b[0m \u001b[1m mean \u001b[0m \u001b[1m std \u001b[0m \u001b[1m naive_se \u001b[0m \u001b[1m mcse \u001b[0m \u001b[1m ess \u001b[0m \u001b[1m rhat \u001b[0m \u001b[1m e\u001b[0m ⋯\n", " \u001b[90m Symbol \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m \u001b[0m ⋯\n", "\n", " α 0.4584 0.0338 0.0011 0.0021 237.0054 0.9998 ⋯\n", " β 0.9516 0.0144 0.0005 0.0009 212.9579 1.0014 ⋯\n", "\u001b[36m 1 column omitted\u001b[0m\n", "\n", "Quantiles\n", " \u001b[1m parameters \u001b[0m \u001b[1m 2.5% \u001b[0m \u001b[1m 25.0% \u001b[0m \u001b[1m 50.0% \u001b[0m \u001b[1m 75.0% \u001b[0m \u001b[1m 97.5% \u001b[0m\n", " \u001b[90m Symbol \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m\n", "\n", " α 0.3922 0.4357 0.4582 0.4806 0.5219\n", " β 0.9212 0.9424 0.9530 0.9614 0.9771\n" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "## Estimation example: first-order, marginal likelihood approach\n", "\n", "# Turing model definition\n", "@model function rbc_kalman(z, m, p_f, cache)\n", " α ~ Uniform(0.2, 0.8)\n", " β ~ Uniform(0.5, 0.99)\n", " p_d = (α = α, β = β)\n", " sol = generate_perturbation(m, p_d, p_f, Val(1); cache)\n", " if !(sol.retcode == :Success)\n", " @addlogprob! -Inf\n", " return\n", " end\n", " @addlogprob! solve(sol, sol.x_ergodic, (0, length(z)); observables = z).logpdf\n", "end\n", "\n", "c = SolverCache(model_rbc, Val(1), p_d)\n", "turing_model = rbc_kalman(fake_z, model_rbc, p_f, c)\n", "n_samples = 1000\n", "n_adapts = 100\n", "δ = 0.65\n", "chain = sample(turing_model, NUTS(n_adapts, δ), n_samples; progress = true)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "┌ Info: Found initial step size\n", "│ ϵ = 0.21250000000000002\n", "└ @ Turing.Inference C:\\Users\\wupei\\.julia\\packages\\Turing\\nfMhU\\src\\inference\\hmc.jl:188\n", "\u001b[32mSampling: 100%|█████████████████████████████████████████| Time: 0:02:26\u001b[39m\n" ] }, { "data": { "text/plain": [ "Chains MCMC chain (1000×35×1 Array{Float64, 3}):\n", "\n", "Iterations = 101:1:1100\n", "Number of chains = 1\n", "Samples per chain = 1000\n", "Wall duration = 182.9 seconds\n", "Compute duration = 182.9 seconds\n", "parameters = α, β, ϵ_draw[1], ϵ_draw[2], ϵ_draw[3], ϵ_draw[4], ϵ_draw[5], ϵ_draw[6], ϵ_draw[7], ϵ_draw[8], ϵ_draw[9], ϵ_draw[10], ϵ_draw[11], ϵ_draw[12], ϵ_draw[13], ϵ_draw[14], ϵ_draw[15], ϵ_draw[16], ϵ_draw[17], ϵ_draw[18], ϵ_draw[19], ϵ_draw[20], ϵ_draw[21]\n", "internals = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size\n", "\n", "Summary Statistics\n", " \u001b[1m parameters \u001b[0m \u001b[1m mean \u001b[0m \u001b[1m std \u001b[0m \u001b[1m naive_se \u001b[0m \u001b[1m mcse \u001b[0m \u001b[1m ess \u001b[0m \u001b[1m rhat \u001b[0m \u001b[1m e\u001b[0m ⋯\n", " \u001b[90m Symbol \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m \u001b[0m ⋯\n", "\n", " α 0.4581 0.0303 0.0010 0.0029 88.5694 1.0073 ⋯\n", " β 0.9530 0.0129 0.0004 0.0010 137.1778 0.9998 ⋯\n", " ϵ_draw[1] -0.1112 0.2148 0.0068 0.0135 254.0827 1.0018 ⋯\n", " ϵ_draw[2] -1.6127 0.3801 0.0120 0.0384 78.9215 1.0020 ⋯\n", " ϵ_draw[3] 0.9781 0.3268 0.0103 0.0275 118.5946 0.9990 ⋯\n", " ϵ_draw[4] 0.0767 0.3007 0.0095 0.0215 191.5072 0.9991 ⋯\n", " ϵ_draw[5] -0.8353 0.3071 0.0097 0.0237 137.3091 1.0009 ⋯\n", " ϵ_draw[6] 0.7644 0.3130 0.0099 0.0212 201.3213 1.0058 ⋯\n", " ϵ_draw[7] -2.1753 0.4250 0.0134 0.0453 55.7146 1.0245 ⋯\n", " ϵ_draw[8] -0.1410 0.3344 0.0106 0.0282 138.6009 1.0063 ⋯\n", " ϵ_draw[9] 0.0578 0.3263 0.0103 0.0282 123.2021 1.0167 ⋯\n", " ϵ_draw[10] 0.3613 0.2978 0.0094 0.0225 150.6162 1.0193 ⋯\n", " ϵ_draw[11] 0.9990 0.3364 0.0106 0.0292 109.8132 1.0318 ⋯\n", " ϵ_draw[12] 0.7090 0.3410 0.0108 0.0311 110.5820 1.0041 ⋯\n", " ϵ_draw[13] 1.3763 0.3569 0.0113 0.0318 121.6980 0.9990 ⋯\n", " ϵ_draw[14] 0.1758 0.3145 0.0099 0.0258 97.5849 1.0293 ⋯\n", " ϵ_draw[15] -1.2367 0.3541 0.0112 0.0360 45.1033 1.0469 ⋯\n", " ϵ_draw[16] -1.4202 0.3882 0.0123 0.0333 132.8966 1.0016 ⋯\n", " ϵ_draw[17] -1.1673 0.3560 0.0113 0.0328 98.0627 1.0097 ⋯\n", " ϵ_draw[18] -0.0066 0.2863 0.0091 0.0192 212.5723 0.9990 ⋯\n", " ϵ_draw[19] 0.6946 0.2987 0.0094 0.0234 134.0651 1.0145 ⋯\n", " ϵ_draw[20] -0.0934 0.3022 0.0096 0.0224 178.9008 1.0169 ⋯\n", " ϵ_draw[21] 0.2455 0.8400 0.0266 0.1138 32.7431 1.0297 ⋯\n", "\u001b[36m 1 column omitted\u001b[0m\n", "\n", "Quantiles\n", " \u001b[1m parameters \u001b[0m \u001b[1m 2.5% \u001b[0m \u001b[1m 25.0% \u001b[0m \u001b[1m 50.0% \u001b[0m \u001b[1m 75.0% \u001b[0m \u001b[1m 97.5% \u001b[0m\n", " \u001b[90m Symbol \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m\n", "\n", " α 0.3986 0.4381 0.4565 0.4796 0.5132\n", " β 0.9272 0.9441 0.9533 0.9627 0.9769\n", " ϵ_draw[1] -0.5821 -0.2391 -0.1085 0.0325 0.2853\n", " ϵ_draw[2] -2.3621 -1.8703 -1.5998 -1.3407 -0.9117\n", " ϵ_draw[3] 0.3892 0.7321 0.9788 1.1906 1.6516\n", " ϵ_draw[4] -0.5363 -0.1241 0.0813 0.2777 0.6609\n", " ϵ_draw[5] -1.4308 -1.0292 -0.8373 -0.6355 -0.2063\n", " ϵ_draw[6] 0.1613 0.5525 0.7625 0.9680 1.3617\n", " ϵ_draw[7] -3.0713 -2.4564 -2.1515 -1.8744 -1.4148\n", " ϵ_draw[8] -0.8838 -0.3437 -0.1170 0.0930 0.4603\n", " ϵ_draw[9] -0.5789 -0.1598 0.0502 0.2642 0.7354\n", " ϵ_draw[10] -0.2095 0.1575 0.3489 0.5470 0.9839\n", " ϵ_draw[11] 0.3521 0.7656 0.9840 1.2140 1.6836\n", " ϵ_draw[12] 0.0783 0.4813 0.6932 0.9084 1.4197\n", " ϵ_draw[13] 0.6989 1.1431 1.3750 1.6057 2.1030\n", " ϵ_draw[14] -0.4388 -0.0387 0.1694 0.3707 0.8161\n", " ϵ_draw[15] -1.9408 -1.4770 -1.2262 -1.0006 -0.5732\n", " ϵ_draw[16] -2.2357 -1.6557 -1.4001 -1.1456 -0.7469\n", " ϵ_draw[17] -1.9231 -1.3939 -1.1531 -0.9090 -0.5508\n", " ϵ_draw[18] -0.5756 -0.1715 -0.0108 0.1741 0.5813\n", " ϵ_draw[19] 0.1203 0.4916 0.6926 0.8879 1.3186\n", " ϵ_draw[20] -0.7047 -0.2895 -0.0822 0.1252 0.4489\n", " ϵ_draw[21] -1.4493 -0.2894 0.2262 0.8156 1.9122\n" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "## Estimation example: first-order, joint likelihood approach\n", "\n", "# Turing model definition\n", "@model function rbc_joint(z, m, p_f, cache, x0 = zeros(m.n_x))\n", " α ~ Uniform(0.2, 0.8)\n", " β ~ Uniform(0.5, 0.99)\n", " p_d = (α = α, β = β)\n", " T = length(z)\n", " ϵ_draw ~ MvNormal(T, 1.0)\n", " ϵ = map(i -> ϵ_draw[((i-1)*m.n_ϵ+1):(i*m.n_ϵ)], 1:T)\n", " # println(p_d)\n", " sol = generate_perturbation(m, p_d, p_f, Val(1); cache)\n", " if !(sol.retcode == :Success)\n", " @addlogprob! -Inf\n", " return\n", " end\n", " @addlogprob! solve(sol, x0, (0, T); noise = ϵ, observables = z).logpdf\n", "end\n", "\n", "c = SolverCache(model_rbc, Val(1), p_d)\n", "turing_model = rbc_joint(fake_z, model_rbc, p_f, c)\n", "n_samples = 1000\n", "n_adapts = 100\n", "δ = 0.65\n", "max_depth = 5 # A lower max_depth will lead to higher autocorrelation of samples, but faster. The time complexity is approximately 2^max_depth\n", "chain = sample(turing_model, NUTS(n_adapts, δ; max_depth), n_samples; progress = true)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "┌ Info: Found initial step size\n", "│ ϵ = 0.3361328125\n", "└ @ Turing.Inference C:\\Users\\wupei\\.julia\\packages\\Turing\\nfMhU\\src\\inference\\hmc.jl:188\n", "\u001b[32mSampling: 100%|█████████████████████████████████████████| Time: 0:04:48\u001b[39m\n" ] }, { "data": { "text/plain": [ "Chains MCMC chain (1000×35×1 Array{Float64, 3}):\n", "\n", "Iterations = 101:1:1100\n", "Number of chains = 1\n", "Samples per chain = 1000\n", "Wall duration = 314.0 seconds\n", "Compute duration = 314.0 seconds\n", "parameters = α, β, ϵ_draw[1], ϵ_draw[2], ϵ_draw[3], ϵ_draw[4], ϵ_draw[5], ϵ_draw[6], ϵ_draw[7], ϵ_draw[8], ϵ_draw[9], ϵ_draw[10], ϵ_draw[11], ϵ_draw[12], ϵ_draw[13], ϵ_draw[14], ϵ_draw[15], ϵ_draw[16], ϵ_draw[17], ϵ_draw[18], ϵ_draw[19], ϵ_draw[20], ϵ_draw[21]\n", "internals = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size\n", "\n", "Summary Statistics\n", " \u001b[1m parameters \u001b[0m \u001b[1m mean \u001b[0m \u001b[1m std \u001b[0m \u001b[1m naive_se \u001b[0m \u001b[1m mcse \u001b[0m \u001b[1m ess \u001b[0m \u001b[1m rhat \u001b[0m \u001b[1m e\u001b[0m ⋯\n", " \u001b[90m Symbol \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m \u001b[0m ⋯\n", "\n", " α 0.4464 0.0365 0.0012 0.0048 21.0310 1.0559 ⋯\n", " β 0.9576 0.0163 0.0005 0.0020 23.2836 1.0623 ⋯\n", " ϵ_draw[1] -0.1132 0.2236 0.0071 0.0129 322.5306 1.0038 ⋯\n", " ϵ_draw[2] -1.6483 0.3299 0.0104 0.0257 135.1192 1.0057 ⋯\n", " ϵ_draw[3] 0.9946 0.3257 0.0103 0.0228 162.2103 0.9990 ⋯\n", " ϵ_draw[4] 0.0842 0.3149 0.0100 0.0204 195.4256 1.0053 ⋯\n", " ϵ_draw[5] -0.8197 0.3311 0.0105 0.0244 169.4104 0.9991 ⋯\n", " ϵ_draw[6] 0.7392 0.3197 0.0101 0.0235 180.6771 0.9997 ⋯\n", " ϵ_draw[7] -2.1800 0.3829 0.0121 0.0353 94.6605 0.9994 ⋯\n", " ϵ_draw[8] -0.1633 0.3227 0.0102 0.0218 209.1768 0.9991 ⋯\n", " ϵ_draw[9] 0.0988 0.3416 0.0108 0.0276 148.3484 0.9994 ⋯\n", " ϵ_draw[10] 0.3421 0.3350 0.0106 0.0224 192.7470 0.9995 ⋯\n", " ϵ_draw[11] 1.0727 0.3705 0.0117 0.0339 124.0563 0.9991 ⋯\n", " ϵ_draw[12] 0.6768 0.3368 0.0106 0.0254 166.8420 1.0030 ⋯\n", " ϵ_draw[13] 1.4199 0.3660 0.0116 0.0297 143.1292 0.9990 ⋯\n", " ϵ_draw[14] 0.1841 0.3071 0.0097 0.0211 189.1216 1.0200 ⋯\n", " ϵ_draw[15] -1.2584 0.3282 0.0104 0.0214 207.2431 1.0041 ⋯\n", " ϵ_draw[16] -1.4749 0.3792 0.0120 0.0325 117.7249 1.0022 ⋯\n", " ϵ_draw[17] -1.1119 0.3349 0.0106 0.0248 155.6362 1.0149 ⋯\n", " ϵ_draw[18] -0.0729 0.3063 0.0097 0.0241 157.5904 1.0055 ⋯\n", " ϵ_draw[19] 0.7628 0.3302 0.0104 0.0266 153.4378 0.9990 ⋯\n", " ϵ_draw[20] -0.1290 0.3373 0.0107 0.0244 181.6521 0.9993 ⋯\n", " ϵ_draw[21] 0.1036 1.0803 0.0342 0.1573 26.9836 0.9996 ⋯\n", "\u001b[36m 1 column omitted\u001b[0m\n", "\n", "Quantiles\n", " \u001b[1m parameters \u001b[0m \u001b[1m 2.5% \u001b[0m \u001b[1m 25.0% \u001b[0m \u001b[1m 50.0% \u001b[0m \u001b[1m 75.0% \u001b[0m \u001b[1m 97.5% \u001b[0m\n", " \u001b[90m Symbol \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m\n", "\n", " α 0.3680 0.4224 0.4504 0.4719 0.5092\n", " β 0.9292 0.9452 0.9559 0.9681 0.9890\n", " ϵ_draw[1] -0.5377 -0.2575 -0.1143 0.0366 0.3093\n", " ϵ_draw[2] -2.3022 -1.8760 -1.6345 -1.4229 -1.0078\n", " ϵ_draw[3] 0.3920 0.7796 0.9840 1.2017 1.6058\n", " ϵ_draw[4] -0.4895 -0.1231 0.0836 0.2986 0.7042\n", " ϵ_draw[5] -1.4837 -1.0168 -0.8432 -0.6017 -0.1574\n", " ϵ_draw[6] 0.1416 0.5062 0.7398 0.9489 1.3727\n", " ϵ_draw[7] -2.9379 -2.4387 -2.1843 -1.9008 -1.4795\n", " ϵ_draw[8] -0.8892 -0.3809 -0.1370 0.0618 0.4209\n", " ϵ_draw[9] -0.5651 -0.1390 0.1001 0.3227 0.7729\n", " ϵ_draw[10] -0.2854 0.1179 0.3234 0.5570 1.0664\n", " ϵ_draw[11] 0.4102 0.8309 1.0583 1.3086 1.8208\n", " ϵ_draw[12] 0.0687 0.4465 0.6619 0.8864 1.3733\n", " ϵ_draw[13] 0.7310 1.1664 1.4166 1.6612 2.1622\n", " ϵ_draw[14] -0.4114 -0.0211 0.2044 0.3867 0.7812\n", " ϵ_draw[15] -1.9786 -1.4656 -1.2507 -1.0293 -0.6627\n", " ϵ_draw[16] -2.2371 -1.7418 -1.4459 -1.2001 -0.8466\n", " ϵ_draw[17] -1.7675 -1.3303 -1.0980 -0.8747 -0.4988\n", " ϵ_draw[18] -0.6794 -0.2647 -0.0651 0.1333 0.5072\n", " ϵ_draw[19] 0.1462 0.5362 0.7524 0.9908 1.3992\n", " ϵ_draw[20] -0.7910 -0.3529 -0.1329 0.0879 0.5805\n", " ϵ_draw[21] -1.9288 -0.7344 0.1410 0.9009 2.3125\n" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "## Estimation example: second-order, joint likelihood approach\n", "\n", "# Turing model definition\n", "@model function rbc_second(z, m, p_f, cache, x0 = zeros(m.n_x))\n", " α ~ Uniform(0.2, 0.8)\n", " β ~ Uniform(0.5, 0.99)\n", " p_d = (α = α, β = β)\n", " T = length(z)\n", " ϵ_draw ~ MvNormal(T, 1.0)\n", " ϵ = map(i -> ϵ_draw[((i-1)*m.n_ϵ+1):(i*m.n_ϵ)], 1:T)\n", " sol = generate_perturbation(m, p_d, p_f, Val(2); cache)\n", " if !(sol.retcode == :Success)\n", " @addlogprob! -Inf\n", " return\n", " end\n", " @addlogprob! solve(sol, x0, (0, T); noise = ϵ, observables = z).logpdf\n", "end\n", "\n", "c = SolverCache(model_rbc, Val(2), p_d)\n", "turing_model = rbc_second(fake_z_second, model_rbc, p_f, c)\n", "n_samples = 1000\n", "n_adapts = 100\n", "δ = 0.65\n", "max_depth = 5\n", "chain = sample(turing_model, NUTS(n_adapts, δ; max_depth), n_samples; progress = true)" ] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.6.4", "language": "julia", "name": "julia-1.6" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.6.4" } }, "nbformat": 4, "nbformat_minor": 4 }