<!-- markdown css tag --><div class="pinggu_markdown">
<div class="pinggu_markdown__html"><blockquote>
<p>DSGE Particle Filtering 粒子滤波<br>
School:Renmin University of China. 中国人民大学 经济学院<br>
Author: Daniel tulips liu</p>
</blockquote>
<h2 id="copyright-c-2&#48;13-2&#48;15-benjamin-born--johannes-pfeifer"><em><strong>Copyright © 2&#48;13-2&#48;15 Benjamin Born + Johannes Pfeifer</strong></em></h2>
<h2 id="main-program---技巧在程序语言上设定为c语言其实是matlab）">main program (技巧在程序语言上设定为C语言其实是MATLAB）</h2>
<p><em><strong>run_filter_and_smoother_AR1_TaRB</strong></em></p>
<pre class=" language-c"><code class="prismlanguage-c"><span class="token operator">%</span> Generate Data from AR1<span class="token operator">-</span>stochastic volatility model<span class="token punctuation">,</span> run MCMC using particle filter on it<span class="token punctuation">.</span> Then use smoother<span class="token punctuation">.</span>
<span class="token operator">%</span>
<span class="token operator">%</span> Notes<span class="token punctuation">:</span>
<span class="token operator">%</span> <span class="token operator">-</span> This code does not <span class="token keyword">do</span> mode<span class="token operator">-</span>finding and computing Hessian at the mode<span class="token punctuation">,</span> because the likelihood from the particle filter is not differentiable<span class="token punctuation">.</span>
<span class="token operator">%</span> Instead<span class="token punctuation">,</span> a diagonal matrix is used<span class="token punctuation">.</span> This implies loss of efficiency but any positive definite matrix will still yield draws from the ergodic distribution
<span class="token operator">%</span> <span class="token punctuation">(</span>see Chib<span class="token operator">/</span><span class="token function">Greenberg</span> <span class="token punctuation">(</span><span class="token number">1995</span><span class="token punctuation">)</span><span class="token punctuation">)</span>
<span class="token operator">%</span> <span class="token operator">-</span> To perform mode<span class="token operator">-</span>finding using the CMA<span class="token operator">-</span>ES algorithm<span class="token punctuation">,</span> uncomment the code in lines <span class="token number">1&#48;1</span> following and
<span class="token operator">%</span> download the cmaes<span class="token punctuation">.</span>m from the provided link<span class="token punctuation">.</span>
<span class="token operator">%</span> <span class="token operator">-</span> Using an insufficient number of particles will yield problems with the acceptance rate<span class="token punctuation">,</span> see Pitt et al<span class="token punctuation">.</span> <span class="token punctuation">(</span><span class="token number">2&#48;12</span><span class="token punctuation">)</span><span class="token punctuation">:</span>
<span class="token operator">%</span> <span class="token string">"On some properties of Markov chain Monte Carlo simulation methods based on the particle filter"</span><span class="token punctuation">,</span> Journal of Econometrics<span class="token punctuation">,</span> <span class="token number">171</span><span class="token punctuation">,</span> <span class="token number">134</span><span class="token operator">-</span><span class="token number">151</span>
<span class="token operator">%</span> <span class="token operator">-</span> In contrast to SMM approaches<span class="token punctuation">,</span> the simulations should not be conducted with fixed random numbers<span class="token punctuation">,</span> unless the number of particles is really large<span class="token punctuation">.</span> See
<span class="token operator">%</span> Thomas Flury and Neil <span class="token function">Shephard</span> <span class="token punctuation">(</span><span class="token number">2&#48;11</span><span class="token punctuation">)</span><span class="token punctuation">:</span> <span class="token string">"Bayesian Inference based only on simulated likelihood"</span><span class="token punctuation">,</span> Econometric Theory<span class="token punctuation">,</span> <span class="token number">27</span><span class="token punctuation">,</span> <span class="token number">933</span>�<span class="token number">956</span>
<span class="token operator">%</span>
<span class="token operator">%</span> <span class="token function">Copyright</span> <span class="token punctuation">(</span>C<span class="token punctuation">)</span> <span class="token number">2&#48;13</span><span class="token operator">-</span><span class="token number">2&#48;15</span> Benjamin Born <span class="token operator">+</span> Johannes Pfeifer
<span class="token operator">%</span>
<span class="token operator">%</span> This is free software<span class="token punctuation">:</span> you can redistribute it and<span class="token operator">/</span>or modify
<span class="token operator">%</span> it under the terms of the GNU General Public License as published by
<span class="token operator">%</span> the Free Software Foundation<span class="token punctuation">,</span> either version <span class="token number">3</span> of the License<span class="token punctuation">,</span> or
<span class="token operator">%</span> <span class="token punctuation">(</span>at your option<span class="token punctuation">)</span> any later version<span class="token punctuation">.</span>
<span class="token operator">%</span>
<span class="token operator">%</span> It is distributed in the hope that it will be useful<span class="token punctuation">,</span>
<span class="token operator">%</span> but WITHOUT ANY WARRANTY<span class="token punctuation">;</span> without even the implied warranty of
<span class="token operator">%</span> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE<span class="token punctuation">.</span>See the
<span class="token operator">%</span> GNU General Public License <span class="token keyword">for</span> more details<span class="token punctuation">.</span>
<span class="token operator">%</span>
<span class="token operator">%</span> See <span class="token operator">&lt;</span>http<span class="token punctuation">:</span><span class="token operator">/</span><span class="token operator">/</span>www<span class="token punctuation">.</span>gnu<span class="token punctuation">.</span>org<span class="token operator">/</span>licenses<span class="token operator">/</span><span class="token operator">&gt;</span><span class="token punctuation">.</span>

close all
clear all
clc

<span class="token function">randn</span><span class="token punctuation">(</span><span class="token string">'state'</span><span class="token punctuation">,</span><span class="token number">1</span><span class="token punctuation">)</span>
<span class="token function">addpath</span><span class="token punctuation">(</span><span class="token string">'Tools'</span><span class="token punctuation">)</span>
<span class="token operator">%</span><span class="token operator">%</span> generate data
periods<span class="token operator">=</span><span class="token number">1&#48;&#48;&#48;</span><span class="token punctuation">;</span>
rho<span class="token operator">=</span><span class="token number">&#48;.95</span><span class="token punctuation">;</span>
sigma_bar<span class="token operator">=</span><span class="token function">log</span><span class="token punctuation">(</span><span class="token number">&#48;.&#48;1</span><span class="token punctuation">)</span><span class="token punctuation">;</span>
rho_sigma<span class="token operator">=</span><span class="token number">&#48;.9</span><span class="token punctuation">;</span>
eta_sigma<span class="token operator">=</span><span class="token number">&#48;.3</span><span class="token punctuation">;</span>

siggma<span class="token operator">=</span><span class="token function">zeros</span><span class="token punctuation">(</span>periods<span class="token punctuation">,</span><span class="token number">1</span><span class="token punctuation">)</span><span class="token punctuation">;</span>
x<span class="token operator">=</span><span class="token function">zeros</span><span class="token punctuation">(</span>periods<span class="token punctuation">,</span><span class="token number">1</span><span class="token punctuation">)</span><span class="token punctuation">;</span>
<span class="token function">siggma</span><span class="token punctuation">(</span><span class="token number">1</span><span class="token punctuation">)</span><span class="token operator">=</span>sigma_bar<span class="token punctuation">;</span>
epsil<span class="token operator">=</span><span class="token function">randn</span><span class="token punctuation">(</span>periods<span class="token punctuation">,</span><span class="token number">2</span><span class="token punctuation">)</span><span class="token punctuation">;</span>
<span class="token keyword">for</span> jj<span class="token operator">=</span><span class="token number">2</span><span class="token punctuation">:</span>periods
<span class="token function">siggma</span><span class="token punctuation">(</span>jj<span class="token punctuation">)</span><span class="token operator">=</span><span class="token punctuation">(</span><span class="token number">1</span><span class="token operator">-</span>rho_sigma<span class="token punctuation">)</span><span class="token operator">*</span>sigma_bar<span class="token operator">+</span>rho_sigma<span class="token operator">*</span><span class="token function">siggma</span><span class="token punctuation">(</span>jj<span class="token number">-1</span><span class="token punctuation">)</span><span class="token operator">+</span>eta_sigma<span class="token operator">*</span><span class="token function">epsil</span><span class="token punctuation">(</span>jj<span class="token punctuation">,</span><span class="token number">1</span><span class="token punctuation">)</span><span class="token punctuation">;</span>
<span class="token function">x</span><span class="token punctuation">(</span>jj<span class="token punctuation">)</span><span class="token operator">=</span>rho<span class="token operator">*</span><span class="token function">x</span><span class="token punctuation">(</span>jj<span class="token number">-1</span><span class="token punctuation">)</span><span class="token operator">+</span><span class="token function">exp</span><span class="token punctuation">(</span><span class="token function">siggma</span><span class="token punctuation">(</span>jj<span class="token punctuation">)</span><span class="token punctuation">)</span><span class="token operator">*</span><span class="token function">epsil</span><span class="token punctuation">(</span>jj<span class="token punctuation">,</span><span class="token number">2</span><span class="token punctuation">)</span><span class="token punctuation">;</span>
end
true_par<span class="token operator">=</span><span class="token punctuation">[</span>rho_sigma<span class="token punctuation">,</span>rho<span class="token punctuation">,</span>eta_sigma<span class="token punctuation">,</span>sigma_bar<span class="token punctuation">]</span><span class="token punctuation">;</span>
save data_simulated x siggma rho sigma_bar rho_sigma eta_sigma true_par


<span class="token operator">%</span><span class="token operator">%</span> settings

num_sim_filter <span class="token operator">=</span> <span class="token number">1&#48;&#48;&#48;</span><span class="token punctuation">;</span> <span class="token operator">%</span> Number of particles <span class="token keyword">for</span> a particle filter
num_sim_smoother <span class="token operator">=</span> <span class="token number">1&#48;&#48;&#48;</span><span class="token punctuation">;</span> <span class="token operator">%</span> Number of particles <span class="token keyword">for</span> a particle filter
dimy<span class="token operator">=</span> <span class="token number">1</span><span class="token punctuation">;</span> <span class="token operator">%</span> Number of variables in the set of <span class="token function">observables</span> <span class="token punctuation">(</span>y<span class="token punctuation">)</span>
burnin <span class="token operator">=</span> <span class="token number">5&#48;&#48;</span><span class="token punctuation">;</span>
MH_draws <span class="token operator">=</span> <span class="token number">5&#48;&#48;&#48;</span><span class="token punctuation">;</span>

<span class="token operator">%</span> <span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span>LOAD DATA <span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span><span class="token operator">*</span>
load data_simulated x rho sigma_bar rho_sigma eta_sigma

observable_series<span class="token operator">=</span>x<span class="token punctuation">;</span>
observable_series<span class="token operator">=</span><span class="token function">observable_series</span><span class="token punctuation">(</span><span class="token operator">~</span><span class="token function">isnan</span><span class="token punctuation">(</span>observable_series<span class="token punctuation">)</span><span class="token punctuation">,</span><span class="token punctuation">:</span><span class="token punctuation">)</span><span class="token punctuation">;</span>

<span class="token operator">%</span> transform to unbounded support
<span class="token punctuation">[</span>rho_sigma_trans<span class="token punctuation">,</span>rho_trans<span class="token punctuation">,</span>eta_sigma_trans<span class="token punctuation">,</span>sigma_bar_trans<span class="token punctuation">]</span><span class="token operator">=</span><span class="token function">par_transform_AR1</span><span class="token punctuation">(</span>rho_sigma<span class="token punctuation">,</span>rho<span class="token punctuation">,</span>eta_sigma<span class="token punctuation">,</span>sigma_bar<span class="token punctuation">)</span><span class="token punctuation">;</span>
npar<span class="token operator">=</span><span class="token number">4</span><span class="token punctuation">;</span>

<span class="token operator">%</span><span class="token operator">%</span> initialize random numbers

<span class="token operator">%</span>initialize matrices
draws<span class="token operator">=</span><span class="token function">zeros</span><span class="token punctuation">(</span>MH_draws<span class="token operator">+</span>burnin<span class="token punctuation">,</span>npar<span class="token punctuation">)</span><span class="token punctuation">;</span>
likelihood<span class="token operator">=</span><span class="token function">zeros</span><span class="token punctuation">(</span>MH_draws<span class="token operator">+</span>burnin<span class="token punctuation">,</span><span class="token number">1</span><span class="token punctuation">)</span><span class="token punctuation">;</span>

<span class="token operator">%</span>create save name depending on clock time
tempdate<span class="token operator">=</span><span class="token function">deblank</span><span class="token punctuation">(</span><span class="token function">num2str</span><span class="token punctuation">(</span><span class="token function">fix</span><span class="token punctuation">(</span>clock<span class="token punctuation">)</span><span class="token punctuation">)</span><span class="token punctuation">)</span><span class="token punctuation">;</span>
tempdate<span class="token operator">=</span><span class="token punctuation">[</span><span class="token function">tempdate</span><span class="token punctuation">(</span><span class="token number">1</span><span class="token punctuation">:</span><span class="token number">4</span><span class="token punctuation">)</span><span class="token punctuation">,</span><span class="token string">'_'</span><span class="token punctuation">,</span><span class="token function">strtrim</span><span class="token punctuation">(</span><span class="token function">tempdate</span><span class="token punctuation">(</span><span class="token number">7</span><span class="token punctuation">:</span><span class="token number">1&#48;</span><span class="token punctuation">)</span><span class="token punctuation">)</span><span class="token punctuation">,</span><span class="token string">'_'</span><span class="token punctuation">,</span><span class="token function">strtrim</span><span class="token punctuation">(</span><span class="token function">tempdate</span><span class="token punctuation">(</span><span class="token number">13</span><span class="token punctuation">:</span><span class="token number">17</span><span class="token punctuation">)</span><span class="token punctuation">)</span><span class="token punctuation">,</span><span class="token string">'_'</span><span class="token punctuation">,</span><span class="token function">strtrim</span><span class="token punctuation">(</span><span class="token function">tempdate</span><span class="token punctuation">(</span><span class="token number">2&#48;</span><span class="token punctuation">:</span><span class="token number">24</span><span class="token punctuation">)</span><span class="token punctuation">)</span><span class="token punctuation">,</span><span class="token string">'_'</span><span class="token punctuation">,</span><span class="token function">strtrim</span><span class="token punctuation">(</span><span class="token function">tempdate</span><span class="token punctuation">(</span><span class="token number">27</span><span class="token punctuation">:</span><span class="token number">3&#48;</span><span class="token punctuation">)</span><span class="token punctuation">)</span><span class="token punctuation">]</span><span class="token punctuation">;</span>
fname<span class="token operator">=</span><span class="token punctuation">[</span><span class="token string">'simulation_'</span><span class="token punctuation">,</span>tempdate<span class="token punctuation">]</span><span class="token punctuation">;</span>


<span class="token operator">%</span><span class="token operator">%</span> MCMC

scale_mh <span class="token operator">=</span> <span class="token number">1e-8</span><span class="token punctuation">;</span> <span class="token operator">%</span>scale factor <span class="token keyword">for</span> proposal density to get acceptance rate of <span class="token number">23</span><span class="token operator">-</span><span class="token number">4&#48;</span> percent
accept<span class="token operator">=</span><span class="token number">&#48;</span><span class="token punctuation">;</span>
new_block_probability<span class="token operator">=</span><span class="token number">&#48;.3</span><span class="token punctuation">;</span>

crit<span class="token operator">=</span><span class="token number">1e-4</span><span class="token punctuation">;</span>
nit<span class="token operator">=</span><span class="token number">1&#48;&#48;</span><span class="token punctuation">;</span>

<span class="token operator">%</span>set startig value<span class="token punctuation">,</span> starts at true values and makes mode<span class="token operator">-</span>finding redundant
<span class="token function">draws</span><span class="token punctuation">(</span><span class="token number">1</span><span class="token punctuation">,</span><span class="token punctuation">:</span><span class="token punctuation">)</span><span class="token operator">=</span><span class="token punctuation">[</span>rho_sigma_trans<span class="token punctuation">,</span>rho_trans<span class="token punctuation">,</span>eta_sigma_trans<span class="token punctuation">,</span>sigma_bar_trans<span class="token punctuation">]</span><span class="token punctuation">;</span>

<span class="token operator">%</span> set jumping matrix <span class="token keyword">for</span> proposal density
<span class="token operator">%</span> uses simple identity matrix instead of Hessian at the mode
inv_Hessian <span class="token operator">=</span> <span class="token number">1e-3</span><span class="token operator">*</span><span class="token function">eye</span><span class="token punctuation">(</span><span class="token number">4</span><span class="token punctuation">)</span><span class="token punctuation">;</span>
Sigma_chol <span class="token operator">=</span> <span class="token function">cholcov</span><span class="token punctuation">(</span>inv_Hessian<span class="token punctuation">)</span>'<span class="token punctuation">;</span>

last_posterior <span class="token operator">=</span> <span class="token function">PF_caller</span><span class="token punctuation">(</span><span class="token function">draws</span><span class="token punctuation">(</span><span class="token number">1</span><span class="token punctuation">,</span><span class="token punctuation">:</span><span class="token punctuation">)</span><span class="token punctuation">,</span>observable_series<span class="token punctuation">,</span>num_sim_filter<span class="token punctuation">,</span>num_sim_smoother<span class="token punctuation">,</span><span class="token number">&#48;</span><span class="token punctuation">)</span><span class="token punctuation">;</span>
<span class="token function">likelihood</span><span class="token punctuation">(</span><span class="token number">1</span><span class="token punctuation">,</span><span class="token number">1</span><span class="token punctuation">)</span><span class="token operator">=</span>last_posterior<span class="token punctuation">;</span>

proposal_draws <span class="token operator">=</span> scale_mh<span class="token operator">*</span>Sigma_chol<span class="token operator">*</span><span class="token function">randn</span><span class="token punctuation">(</span><span class="token number">4</span><span class="token punctuation">,</span>MH_draws<span class="token operator">+</span>burnin<span class="token punctuation">)</span><span class="token punctuation">;</span>
blocked_draws_counter<span class="token operator">=</span><span class="token number">&#48;</span><span class="token punctuation">;</span>
accepted_draws_counter<span class="token operator">=</span><span class="token number">&#48;</span><span class="token punctuation">;</span>
<span class="token keyword">for</span> draw_iter<span class="token operator">=</span><span class="token number">2</span><span class="token punctuation">:</span>burnin<span class="token operator">+</span>MH_draws
<span class="token operator">%</span><span class="token operator">%</span> randomize indices <span class="token keyword">for</span> blocking in this iteration
indices<span class="token operator">=</span><span class="token function">randperm</span><span class="token punctuation">(</span>npar<span class="token punctuation">)</span>'<span class="token punctuation">;</span>
blocks<span class="token operator">=</span><span class="token punctuation">[</span><span class="token number">1</span><span class="token punctuation">;</span> <span class="token punctuation">(</span><span class="token number">1</span><span class="token operator">+</span><span class="token function">cumsum</span><span class="token punctuation">(</span><span class="token punctuation">(</span><span class="token function">rand</span><span class="token punctuation">(</span><span class="token function">length</span><span class="token punctuation">(</span>indices<span class="token punctuation">)</span><span class="token operator">-</span><span class="token number">1</span><span class="token punctuation">,</span><span class="token number">1</span><span class="token punctuation">)</span><span class="token operator">&gt;</span><span class="token punctuation">(</span><span class="token number">1</span><span class="token operator">-</span>new_block_probability<span class="token punctuation">)</span><span class="token punctuation">)</span><span class="token punctuation">)</span><span class="token punctuation">)</span><span class="token punctuation">]</span><span class="token punctuation">;</span>
nblocks<span class="token operator">=</span><span class="token function">blocks</span><span class="token punctuation">(</span>end<span class="token punctuation">,</span><span class="token number">1</span><span class="token punctuation">)</span><span class="token punctuation">;</span> <span class="token operator">%</span>get number of blocks this iteration
current_draw<span class="token operator">=</span><span class="token function">draws</span><span class="token punctuation">(</span>draw_iter<span class="token number">-1</span><span class="token punctuation">,</span><span class="token punctuation">:</span><span class="token punctuation">)</span>'<span class="token punctuation">;</span> <span class="token operator">%</span>get starting point <span class="token keyword">for</span> current draw <span class="token keyword">for</span> updating
<span class="token keyword">for</span> block_iter<span class="token operator">=</span><span class="token number">1</span><span class="token punctuation">:</span>nblocks
blocked_draws_counter<span class="token operator">=</span>blocked_draws_counter<span class="token operator">+</span><span class="token number">1</span><span class="token punctuation">;</span>
nxopt<span class="token operator">=</span><span class="token function">length</span><span class="token punctuation">(</span><span class="token function">indices</span><span class="token punctuation">(</span>blocks<span class="token operator">==</span>block_iter<span class="token punctuation">,</span><span class="token number">1</span><span class="token punctuation">)</span><span class="token punctuation">)</span><span class="token punctuation">;</span> <span class="token operator">%</span>get size of current block
par_start_current_block<span class="token operator">=</span><span class="token function">current_draw</span><span class="token punctuation">(</span><span class="token function">indices</span><span class="token punctuation">(</span>blocks<span class="token operator">==</span>block_iter<span class="token punctuation">,</span><span class="token number">1</span><span class="token punctuation">)</span><span class="token punctuation">)</span><span class="token punctuation">;</span>
<span class="token punctuation">[</span>fval<span class="token punctuation">,</span>xopt_current_block<span class="token punctuation">,</span><span class="token operator">~</span><span class="token punctuation">,</span>hessian_mat<span class="token punctuation">,</span><span class="token operator">~</span><span class="token punctuation">,</span><span class="token operator">~</span><span class="token punctuation">,</span>exitflag<span class="token punctuation">]</span> <span class="token operator">=</span> <span class="token punctuation">.</span><span class="token punctuation">.</span><span class="token punctuation">.</span>
<span class="token function">csminwel</span><span class="token punctuation">(</span><span class="token string">'PF_caller_TaRB'</span><span class="token punctuation">,</span>par_start_current_block<span class="token punctuation">,</span><span class="token number">1e-5</span><span class="token operator">*</span><span class="token function">eye</span><span class="token punctuation">(</span>nxopt<span class="token punctuation">)</span><span class="token punctuation">,</span><span class="token punctuation">[</span><span class="token punctuation">]</span><span class="token punctuation">,</span>crit<span class="token punctuation">,</span>nit<span class="token punctuation">,</span><span class="token punctuation">.</span><span class="token punctuation">.</span><span class="token punctuation">.</span>
current_draw<span class="token punctuation">,</span><span class="token function">indices</span><span class="token punctuation">(</span>blocks<span class="token operator">==</span>block_iter<span class="token punctuation">,</span><span class="token number">1</span><span class="token punctuation">)</span><span class="token punctuation">,</span>observable_series<span class="token punctuation">,</span>num_sim_filter<span class="token punctuation">,</span>num_sim_smoother<span class="token punctuation">,</span><span class="token number">&#48;</span><span class="token punctuation">)</span><span class="token punctuation">;</span> <span class="token operator">%</span>inputs <span class="token keyword">for</span> ob jective

<span class="token keyword">if</span> <span class="token function">any</span><span class="token punctuation">(</span><span class="token function">any</span><span class="token punctuation">(</span><span class="token function">isnan</span><span class="token punctuation">(</span>hessian_mat<span class="token punctuation">)</span><span class="token punctuation">)</span><span class="token punctuation">)</span> <span class="token operator">||</span> <span class="token function">any</span><span class="token punctuation">(</span><span class="token function">any</span><span class="token punctuation">(</span><span class="token function">isinf</span><span class="token punctuation">(</span>hessian_mat<span class="token punctuation">)</span><span class="token punctuation">)</span><span class="token punctuation">)</span>
inverse_hessian_mat<span class="token operator">=</span><span class="token function">eye</span><span class="token punctuation">(</span>nxopt<span class="token punctuation">)</span><span class="token operator">*</span><span class="token number">1e-4</span><span class="token punctuation">;</span> <span class="token operator">%</span>use diagonal
<span class="token keyword">else</span>
inverse_hessian_mat<span class="token operator">=</span><span class="token function">inv</span><span class="token punctuation">(</span>hessian_mat<span class="token punctuation">)</span><span class="token punctuation">;</span> <span class="token operator">%</span>get inverse Hessian
<span class="token keyword">if</span> <span class="token function">any</span><span class="token punctuation">(</span><span class="token function">any</span><span class="token punctuation">(</span><span class="token punctuation">(</span><span class="token function">isnan</span><span class="token punctuation">(</span>inverse_hessian_mat<span class="token punctuation">)</span><span class="token punctuation">)</span><span class="token punctuation">)</span><span class="token punctuation">)</span> <span class="token operator">||</span> <span class="token function">any</span><span class="token punctuation">(</span><span class="token function">any</span><span class="token punctuation">(</span><span class="token punctuation">(</span><span class="token function">isinf</span><span class="token punctuation">(</span>inverse_hessian_mat<span class="token punctuation">)</span><span class="token punctuation">)</span><span class="token punctuation">)</span><span class="token punctuation">)</span>
inverse_hessian_mat<span class="token operator">=</span><span class="token function">eye</span><span class="token punctuation">(</span>nxopt<span class="token punctuation">)</span><span class="token operator">*</span><span class="token number">1e-4</span><span class="token punctuation">;</span> <span class="token operator">%</span>use diagonal
end
end
<span class="token punctuation">[</span>proposal_covariance_Cholesky_decomposition_upper<span class="token punctuation">,</span>negeigenvalues<span class="token punctuation">]</span><span class="token operator">=</span><span class="token function">chol</span><span class="token punctuation">(</span>inverse_hessian_mat<span class="token punctuation">)</span><span class="token punctuation">;</span>
<span class="token operator">%</span><span class="token keyword">if</span> not positive definite<span class="token punctuation">,</span> use generalized Cholesky of Eskow<span class="token operator">/</span>Schnabel
<span class="token keyword">if</span> negeigenvalues<span class="token operator">~</span><span class="token operator">=</span><span class="token number">&#48;</span>
proposal_covariance_Cholesky_decomposition_upper<span class="token operator">=</span><span class="token function">chol_SE</span><span class="token punctuation">(</span>inverse_hessian_mat<span class="token punctuation">,</span><span class="token number">&#48;</span><span class="token punctuation">)</span><span class="token punctuation">;</span>
end
proposal_covariance_Cholesky_decomposition_upper<span class="token operator">=</span>proposal_covariance_Cholesky_decomposition_upper<span class="token operator">*</span>scale_mh<span class="token punctuation">;</span>
<span class="token operator">%</span>get proposal draw

tempdraw<span class="token operator">=</span>current_draw<span class="token punctuation">;</span>
r <span class="token operator">=</span> <span class="token function">randn</span><span class="token punctuation">(</span><span class="token number">1</span><span class="token punctuation">,</span>nxopt<span class="token punctuation">)</span> <span class="token operator">*</span> proposal_covariance_Cholesky_decomposition_upper<span class="token punctuation">;</span>
proposed_par<span class="token operator">=</span>xopt_current_block<span class="token operator">+</span>r'<span class="token punctuation">;</span>
<span class="token function">tempdraw</span><span class="token punctuation">(</span><span class="token function">indices</span><span class="token punctuation">(</span>blocks<span class="token operator">==</span>block_iter<span class="token punctuation">,</span><span class="token number">1</span><span class="token punctuation">)</span><span class="token punctuation">)</span><span class="token operator">=</span>proposed_par<span class="token punctuation">;</span>

<span class="token operator">%</span> check whether draw is valid and compute posterior
try
logpost <span class="token operator">=</span> <span class="token operator">-</span> <span class="token punctuation">(</span><span class="token function">PF_caller_TaRB</span><span class="token punctuation">(</span><span class="token function">tempdraw</span><span class="token punctuation">(</span><span class="token function">indices</span><span class="token punctuation">(</span>blocks<span class="token operator">==</span>block_iter<span class="token punctuation">,</span><span class="token number">1</span><span class="token punctuation">)</span><span class="token punctuation">)</span><span class="token punctuation">,</span>current_draw<span class="token punctuation">,</span><span class="token function">indices</span><span class="token punctuation">(</span>blocks<span class="token operator">==</span>block_iter<span class="token punctuation">,</span><span class="token number">1</span><span class="token punctuation">)</span><span class="token punctuation">,</span>observable_series<span class="token punctuation">,</span>num_sim_filter<span class="token punctuation">,</span>num_sim_smoother<span class="token punctuation">,</span><span class="token number">&#48;</span><span class="token punctuation">)</span><span class="token punctuation">)</span><span class="token punctuation">;</span> <span class="token operator">%</span>inputs <span class="token keyword">for</span> ob jective
catch
logpost <span class="token operator">=</span> <span class="token operator">-</span>inf<span class="token punctuation">;</span>
end

<span class="token keyword">if</span> <span class="token punctuation">(</span>logpost <span class="token operator">&gt;</span> <span class="token operator">-</span>inf<span class="token punctuation">)</span>
<span class="token operator">%</span>get ratio of proposal densities<span class="token punctuation">,</span> required because proposal depends
<span class="token operator">%</span>on current mode via Hessian and is thus not symmetric anymore
proposal_density_proposed_move_forward<span class="token operator">=</span><span class="token function">multivariate_normal_pdf</span><span class="token punctuation">(</span>proposed_par<span class="token string">',xopt_current_block'</span><span class="token punctuation">,</span>proposal_covariance_Cholesky_decomposition_upper<span class="token punctuation">,</span>nxopt<span class="token punctuation">)</span><span class="token punctuation">;</span>
proposal_density_proposed_move_backward<span class="token operator">=</span><span class="token function">multivariate_normal_pdf</span><span class="token punctuation">(</span>par_start_current_block<span class="token string">',xopt_current_block'</span><span class="token punctuation">,</span>proposal_covariance_Cholesky_decomposition_upper<span class="token punctuation">,</span>nxopt<span class="token punctuation">)</span><span class="token punctuation">;</span>
accprob<span class="token operator">=</span>logpost<span class="token operator">-</span>last_posterior<span class="token operator">+</span> <span class="token function">log</span><span class="token punctuation">(</span>proposal_density_proposed_move_backward<span class="token punctuation">)</span><span class="token operator">-</span><span class="token function">log</span><span class="token punctuation">(</span>proposal_density_proposed_move_forward<span class="token punctuation">)</span><span class="token punctuation">;</span> <span class="token operator">%</span><span class="token function">Formula</span> <span class="token punctuation">(</span><span class="token number">6</span><span class="token punctuation">)</span><span class="token punctuation">,</span> Chib<span class="token operator">/</span>Ramamurthy
<span class="token function">exp</span><span class="token punctuation">(</span>accprob<span class="token punctuation">)</span>
<span class="token keyword">if</span><span class="token punctuation">(</span><span class="token function">log</span><span class="token punctuation">(</span>rand<span class="token punctuation">)</span> <span class="token operator">&lt;</span> accprob<span class="token punctuation">)</span>
<span class="token function">current_draw</span><span class="token punctuation">(</span><span class="token function">indices</span><span class="token punctuation">(</span>blocks<span class="token operator">==</span>block_iter<span class="token punctuation">,</span><span class="token number">1</span><span class="token punctuation">)</span><span class="token punctuation">)</span><span class="token operator">=</span>proposed_par<span class="token punctuation">;</span>
last_posterior<span class="token operator">=</span>logpost<span class="token punctuation">;</span>
accepted_draws_counter <span class="token operator">=</span>accepted_draws_counter <span class="token operator">+</span><span class="token number">1</span><span class="token punctuation">;</span>
<span class="token keyword">else</span> <span class="token operator">%</span>no updating
<span class="token operator">%</span><span class="token keyword">do</span> nothing<span class="token punctuation">,</span> keep old value
end
end
end
accepted<span class="token operator">=</span>accepted_draws_counter<span class="token operator">/</span>blocked_draws_counter<span class="token punctuation">;</span>
<span class="token function">draws</span><span class="token punctuation">(</span>draw_iter<span class="token punctuation">,</span><span class="token punctuation">:</span><span class="token punctuation">)</span> <span class="token operator">=</span> current_draw<span class="token punctuation">;</span>
logpost <span class="token operator">=</span> last_posterior<span class="token punctuation">;</span> <span class="token operator">%</span>make sure not a temporary draw is returned<span class="token punctuation">;</span>
<span class="token keyword">if</span> <span class="token function">mod</span><span class="token punctuation">(</span>draw_iter<span class="token punctuation">,</span><span class="token number">1&#48;</span><span class="token punctuation">)</span><span class="token operator">==</span><span class="token number">&#48;</span>
<span class="token function">disp</span><span class="token punctuation">(</span><span class="token punctuation">[</span><span class="token string">'Acceptance Rate:'</span> <span class="token function">num2str</span><span class="token punctuation">(</span>accepted<span class="token punctuation">)</span><span class="token punctuation">]</span><span class="token punctuation">)</span>
<span class="token function">eval</span><span class="token punctuation">(</span><span class="token punctuation">[</span><span class="token string">'save '</span><span class="token punctuation">,</span>fname<span class="token punctuation">]</span><span class="token punctuation">)</span><span class="token punctuation">;</span>
end
end
<span class="token operator">%</span> retransform parameters and plot them
<span class="token punctuation">[</span>rho_sigma_vec<span class="token punctuation">,</span>rho_1_vec<span class="token punctuation">,</span>eta_sigma_vec<span class="token punctuation">,</span>sigma_bar_vec<span class="token punctuation">]</span><span class="token operator">=</span><span class="token function">par_retransform_AR1</span><span class="token punctuation">(</span><span class="token function">draws</span><span class="token punctuation">(</span>burnin<span class="token punctuation">:</span>ii<span class="token punctuation">,</span><span class="token punctuation">:</span><span class="token punctuation">)</span><span class="token punctuation">)</span><span class="token punctuation">;</span>
figure
<span class="token function">subplot</span><span class="token punctuation">(</span><span class="token number">4</span><span class="token punctuation">,</span><span class="token number">1</span><span class="token punctuation">,</span><span class="token number">1</span><span class="token punctuation">)</span>
<span class="token function">plot</span><span class="token punctuation">(</span>rho_sigma_vec<span class="token punctuation">)</span>
<span class="token function">title</span><span class="token punctuation">(</span><span class="token string">'\rho_\sigma'</span><span class="token punctuation">)</span>
axis tight
<span class="token function">subplot</span><span class="token punctuation">(</span><span class="token number">4</span><span class="token punctuation">,</span><span class="token number">1</span><span class="token punctuation">,</span><span class="token number">2</span><span class="token punctuation">)</span>
<span class="token function">plot</span><span class="token punctuation">(</span>rho_1_vec<span class="token punctuation">)</span>
<span class="token function">title</span><span class="token punctuation">(</span><span class="token string">'\rho'</span><span class="token punctuation">)</span>
axis tight
<span class="token function">subplot</span><span class="token punctuation">(</span><span class="token number">4</span><span class="token punctuation">,</span><span class="token number">1</span><span class="token punctuation">,</span><span class="token number">3</span><span class="token punctuation">)</span>
<span class="token function">plot</span><span class="token punctuation">(</span>eta_sigma_vec<span class="token punctuation">)</span>
<span class="token function">title</span><span class="token punctuation">(</span><span class="token string">'\eta_\sigma'</span><span class="token punctuation">)</span>
axis tight
<span class="token function">subplot</span><span class="token punctuation">(</span><span class="token number">4</span><span class="token punctuation">,</span><span class="token number">1</span><span class="token punctuation">,</span><span class="token number">4</span><span class="token punctuation">)</span>
<span class="token function">plot</span><span class="token punctuation">(</span>sigma_bar_vec<span class="token punctuation">)</span>
<span class="token function">title</span><span class="token punctuation">(</span><span class="token string">'\bar \sigma'</span><span class="token punctuation">)</span>
axis tight

<span class="token function">eval</span><span class="token punctuation">(</span><span class="token punctuation">[</span><span class="token string">'save '</span><span class="token punctuation">,</span>fname<span class="token punctuation">]</span><span class="token punctuation">)</span><span class="token punctuation">;</span>

<span class="token operator">%</span><span class="token operator">%</span> Run Smoother
num_sim_smoother <span class="token operator">=</span> <span class="token number">2&#48;&#48;</span><span class="token punctuation">;</span>

<span class="token operator">%</span> comment the following four line out to obtain smoothed estimates at the
<span class="token operator">%</span> true values

rho_sigma_trans<span class="token operator">=</span><span class="token function">mean</span><span class="token punctuation">(</span><span class="token function">draws</span><span class="token punctuation">(</span>burnin<span class="token punctuation">:</span>ii<span class="token punctuation">,</span><span class="token number">1</span><span class="token punctuation">)</span><span class="token punctuation">)</span><span class="token punctuation">;</span>
rho_trans<span class="token operator">=</span><span class="token function">mean</span><span class="token punctuation">(</span><span class="token function">draws</span><span class="token punctuation">(</span>burnin<span class="token punctuation">:</span>ii<span class="token punctuation">,</span><span class="token number">2</span><span class="token punctuation">)</span><span class="token punctuation">)</span><span class="token punctuation">;</span>
eta_sigma_trans<span class="token operator">=</span><span class="token function">mean</span><span class="token punctuation">(</span><span class="token function">draws</span><span class="token punctuation">(</span>burnin<span class="token punctuation">:</span>ii<span class="token punctuation">,</span><span class="token number">3</span><span class="token punctuation">)</span><span class="token punctuation">)</span><span class="token punctuation">;</span>
sigma_bar_trans<span class="token operator">=</span><span class="token function">mean</span><span class="token punctuation">(</span><span class="token function">draws</span><span class="token punctuation">(</span>burnin<span class="token punctuation">:</span>ii<span class="token punctuation">,</span><span class="token number">4</span><span class="token punctuation">)</span><span class="token punctuation">)</span><span class="token punctuation">;</span>

<span class="token operator">%</span> make sure to use parameters as they are the input arguments<span class="token operator">!</span>
<span class="token function">PF_caller</span><span class="token punctuation">(</span><span class="token punctuation">[</span>rho_sigma_trans<span class="token punctuation">,</span>rho_trans<span class="token punctuation">,</span>eta_sigma_trans<span class="token punctuation">,</span>sigma_bar_trans<span class="token punctuation">]</span><span class="token punctuation">,</span>observable_series<span class="token punctuation">,</span>num_sim_filter<span class="token punctuation">,</span>num_sim_smoother<span class="token punctuation">,</span><span class="token number">1</span><span class="token punctuation">,</span><span class="token string">'Stoch_vol_AR_smoother'</span><span class="token punctuation">)</span><span class="token punctuation">;</span>

<span class="token operator">%</span> load smoother results
load Stoch_vol_AR_smoother vola_resids level_resids x_smoother_store xhat_smoother x_resampled_store

<span class="token operator">%</span>plot them
figure
<span class="token function">subplot</span><span class="token punctuation">(</span><span class="token number">3</span><span class="token punctuation">,</span><span class="token number">1</span><span class="token punctuation">,</span><span class="token number">1</span><span class="token punctuation">)</span>
<span class="token function">plot</span><span class="token punctuation">(</span><span class="token number">1</span><span class="token punctuation">:</span>periods<span class="token punctuation">,</span>xhat_smoother<span class="token punctuation">,</span><span class="token string">'r--'</span><span class="token punctuation">,</span><span class="token string">'LineWidth'</span><span class="token punctuation">,</span><span class="token number">&#48;.5</span><span class="token punctuation">)</span>
<span class="token function">title</span><span class="token punctuation">(</span><span class="token string">'Volatility'</span><span class="token punctuation">)</span>
hold on
<span class="token function">plot</span><span class="token punctuation">(</span><span class="token number">1</span><span class="token punctuation">:</span>periods<span class="token punctuation">,</span>siggma<span class="token punctuation">,</span><span class="token string">'b-'</span><span class="token punctuation">,</span><span class="token string">'LineWidth'</span><span class="token punctuation">,</span><span class="token number">&#48;.5</span><span class="token punctuation">)</span>
<span class="token function">plot</span><span class="token punctuation">(</span><span class="token number">1</span><span class="token punctuation">:</span>periods<span class="token punctuation">,</span><span class="token function">squeeze</span><span class="token punctuation">(</span><span class="token function">median</span><span class="token punctuation">(</span>x_resampled_store<span class="token punctuation">,</span><span class="token number">2</span><span class="token punctuation">)</span><span class="token punctuation">)</span><span class="token punctuation">,</span><span class="token string">'c-.'</span><span class="token punctuation">,</span><span class="token string">'LineWidth'</span><span class="token punctuation">,</span><span class="token number">&#48;.5</span><span class="token punctuation">)</span>
<span class="token function">legend</span><span class="token punctuation">(</span><span class="token string">'Smoothed'</span><span class="token punctuation">,</span><span class="token string">'True'</span><span class="token punctuation">,</span><span class="token string">'Filtered'</span><span class="token punctuation">,</span><span class="token string">'Location'</span><span class="token punctuation">,</span><span class="token string">'NorthWest'</span><span class="token punctuation">)</span>
<span class="token function">fprintf</span><span class="token punctuation">(</span><span class="token string">'Correlation States: %4.3f\n'</span><span class="token punctuation">,</span><span class="token function">corr</span><span class="token punctuation">(</span>xhat_smoother'<span class="token punctuation">,</span>siggma<span class="token punctuation">)</span><span class="token punctuation">)</span>

<span class="token function">subplot</span><span class="token punctuation">(</span><span class="token number">3</span><span class="token punctuation">,</span><span class="token number">1</span><span class="token punctuation">,</span><span class="token number">2</span><span class="token punctuation">)</span>
<span class="token function">plot</span><span class="token punctuation">(</span><span class="token number">1</span><span class="token punctuation">:</span>periods<span class="token punctuation">,</span><span class="token function">median</span><span class="token punctuation">(</span>vola_resids<span class="token punctuation">,</span><span class="token number">2</span><span class="token punctuation">)</span><span class="token operator">/</span>eta_sigma<span class="token punctuation">,</span><span class="token string">'r--'</span><span class="token punctuation">,</span><span class="token string">'LineWidth'</span><span class="token punctuation">,</span><span class="token number">&#48;.5</span><span class="token punctuation">)</span>
hold on
<span class="token function">plot</span><span class="token punctuation">(</span><span class="token number">1</span><span class="token punctuation">:</span>periods<span class="token punctuation">,</span>eta_sigma<span class="token operator">*</span><span class="token function">epsil</span><span class="token punctuation">(</span><span class="token punctuation">:</span><span class="token punctuation">,</span><span class="token number">1</span><span class="token punctuation">)</span><span class="token punctuation">,</span><span class="token string">'b-'</span><span class="token punctuation">,</span><span class="token string">'LineWidth'</span><span class="token punctuation">,</span><span class="token number">&#48;.5</span><span class="token punctuation">)</span>
<span class="token function">title</span><span class="token punctuation">(</span><span class="token string">'Vola Residuals'</span><span class="token punctuation">)</span>
<span class="token function">fprintf</span><span class="token punctuation">(</span><span class="token string">'Correlation State Resids: %4.3f\n'</span><span class="token punctuation">,</span><span class="token function">corr</span><span class="token punctuation">(</span><span class="token function">median</span><span class="token punctuation">(</span>vola_resids<span class="token punctuation">,</span><span class="token number">2</span><span class="token punctuation">)</span><span class="token punctuation">,</span>eta_sigma<span class="token operator">*</span><span class="token function">epsil</span><span class="token punctuation">(</span><span class="token punctuation">:</span><span class="token punctuation">,</span><span class="token number">1</span><span class="token punctuation">)</span><span class="token punctuation">)</span><span class="token punctuation">)</span>
<span class="token function">fprintf</span><span class="token punctuation">(</span><span class="token string">'Std State Resids: %4.3f\n'</span><span class="token punctuation">,</span><span class="token function">std</span><span class="token punctuation">(</span><span class="token function">vola_resids</span><span class="token punctuation">(</span><span class="token punctuation">:</span><span class="token punctuation">)</span><span class="token operator">/</span>eta_sigma<span class="token punctuation">)</span><span class="token punctuation">)</span>

<span class="token function">subplot</span><span class="token punctuation">(</span><span class="token number">3</span><span class="token punctuation">,</span><span class="token number">1</span><span class="token punctuation">,</span><span class="token number">3</span><span class="token punctuation">)</span>
<span class="token function">plot</span><span class="token punctuation">(</span><span class="token number">1</span><span class="token punctuation">:</span>periods<span class="token punctuation">,</span><span class="token function">median</span><span class="token punctuation">(</span>level_resids<span class="token punctuation">,</span><span class="token number">2</span><span class="token punctuation">)</span><span class="token punctuation">,</span><span class="token string">'r--'</span><span class="token punctuation">,</span><span class="token string">'LineWidth'</span><span class="token punctuation">,</span><span class="token number">&#48;.5</span><span class="token punctuation">)</span>
hold on
<span class="token function">plot</span><span class="token punctuation">(</span><span class="token number">1</span><span class="token punctuation">:</span>periods<span class="token punctuation">,</span><span class="token function">epsil</span><span class="token punctuation">(</span><span class="token punctuation">:</span><span class="token punctuation">,</span><span class="token number">2</span><span class="token punctuation">)</span><span class="token punctuation">,</span><span class="token string">'b-'</span><span class="token punctuation">,</span><span class="token string">'LineWidth'</span><span class="token punctuation">,</span><span class="token number">&#48;.5</span><span class="token punctuation">)</span>
<span class="token function">title</span><span class="token punctuation">(</span><span class="token string">'Level Residuals'</span><span class="token punctuation">)</span>
<span class="token function">fprintf</span><span class="token punctuation">(</span><span class="token string">'Correlation Level Resids: %4.3f\n'</span><span class="token punctuation">,</span><span class="token function">corr</span><span class="token punctuation">(</span><span class="token function">median</span><span class="token punctuation">(</span>level_resids<span class="token punctuation">,</span><span class="token number">2</span><span class="token punctuation">)</span><span class="token punctuation">,</span><span class="token function">epsil</span><span class="token punctuation">(</span><span class="token punctuation">:</span><span class="token punctuation">,</span><span class="token number">2</span><span class="token punctuation">)</span><span class="token punctuation">)</span><span class="token punctuation">)</span>
<span class="token function">fprintf</span><span class="token punctuation">(</span><span class="token string">'Std Level Resids: %4.3f\n'</span><span class="token punctuation">,</span><span class="token function">std</span><span class="token punctuation">(</span><span class="token function">level_resids</span><span class="token punctuation">(</span><span class="token punctuation">:</span><span class="token punctuation">)</span><span class="token punctuation">)</span><span class="token punctuation">)</span>

</code></pre>
<h2 id="压缩文件在第一楼收取一点论坛币">压缩文件在第一楼收取一点论坛币</h2>
<p>看一楼</p>
<h2 id="参考资料">参考资料</h2>
<p>[1] dynare team , Pfeifer.</p>
</div>
</div>