Skip to content
GitLab
Menu
Projects
Groups
Snippets
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in / Register
Toggle navigation
Menu
Open sidebar
Mark Song
particles
Commits
4921000f
Commit
4921000f
authored
May 18, 2017
by
Stéphane Adjemian
Browse files
Fixed indentation.
parent
823d9474
Changes
26
Hide whitespace changes
Inline
Side-by-side
src/auxiliary_initialization.m
View file @
4921000f
...
...
@@ -87,7 +87,7 @@ yhat = bsxfun(@minus,StateVectors,state_variables_steady_state);
% yhat_ = bsxfun(@minus,StateVectors_,state_variables_steady_state);
% [tmp, tmp_] = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,ThreadsOptions.local_state_space_iteration_2);
%else
tmp
=
local_state_space_iteration_2
(
yhat
,
zeros
(
number_of_structural_innovations
,
number_of_particles
),
ghx
,
ghu
,
constant
,
ghxx
,
ghuu
,
ghxu
,
ThreadsOptions
.
local_state_space_iteration_2
);
tmp
=
local_state_space_iteration_2
(
yhat
,
zeros
(
number_of_structural_innovations
,
number_of_particles
),
ghx
,
ghu
,
constant
,
ghxx
,
ghuu
,
ghxu
,
ThreadsOptions
.
local_state_space_iteration_2
);
%end
PredictedObservedMean
=
weights
*
(
tmp
(
mf1
,:)
'
);
PredictionError
=
bsxfun
(
@
minus
,
Y
(:,
t
),
tmp
(
mf1
,:));
...
...
src/auxiliary_particle_filter.m
View file @
4921000f
function
[
LIK
,
lik
]
=
auxiliary_particle_filter
(
ReducedForm
,
Y
,
start
,
ParticleOptions
,
ThreadsOptions
)
% Evaluates the likelihood of a nonlinear model with the auxiliary particle filter
% Evaluates the likelihood of a nonlinear model with the auxiliary particle filter
% allowing eventually resampling.
%
% Copyright (C) 2011-2015 Dynare Team
...
...
@@ -91,11 +91,11 @@ end
% [nodes,nodes_weights,nodes_weights_c] = unscented_sigma_points(number_of_structural_innovations,ParticleOptions);
% else
% error('Estimation: This approximation for the proposal is not implemented or unknown!')
% end
% end
% nodes = (Q_lower_triangular_cholesky*(nodes'))' ;
nodes
=
zeros
(
1
,
number_of_structural_innovations
)
;
nodes_weights
=
ones
(
number_of_structural_innovations
,
1
)
;
nodes_weights
=
ones
(
number_of_structural_innovations
,
1
)
;
for
t
=
1
:
sample_size
yhat
=
bsxfun
(
@
minus
,
StateVectors
,
state_variables_steady_state
);
...
...
@@ -126,7 +126,7 @@ for t=1:sample_size
yhat_
=
yhat_
(:,
indx
)
;
end
yhat
=
yhat
(:,
indx
)
;
weights_stage_1
=
weights
(
indx
)
.
/
tau_tilde
(
indx
)
;
weights_stage_1
=
weights
(
indx
)
.
/
tau_tilde
(
indx
)
;
epsilon
=
Q_lower_triangular_cholesky
*
randn
(
number_of_structural_innovations
,
number_of_particles
);
if
pruning
[
tmp
,
tmp_
]
=
local_state_space_iteration_2
(
yhat
,
epsilon
,
ghx
,
ghu
,
constant
,
ghxx
,
ghuu
,
ghxu
,
yhat_
,
steadystate
,
ThreadsOptions
.
local_state_space_iteration_2
);
...
...
@@ -137,7 +137,7 @@ for t=1:sample_size
StateVectors
=
tmp
(
mf0
,:);
PredictionError
=
bsxfun
(
@
minus
,
Y
(:,
t
),
tmp
(
mf1
,:));
weights_stage_2
=
weights_stage_1
.*
(
exp
(
-.
5
*
(
const_lik
+
sum
(
PredictionError
.*
(
H
\
PredictionError
),
1
)))
+
1e-99
)
;
lik
(
t
)
=
log
(
mean
(
weights_stage_2
))
;
lik
(
t
)
=
log
(
mean
(
weights_stage_2
))
;
weights
=
weights_stage_2
/
sum
(
weights_stage_2
);
if
(
ParticleOptions
.
resampling
.
status
.
generic
&&
neff
(
weights
)
<
ParticleOptions
.
resampling
.
threshold
*
sample_size
)
||
ParticleOptions
.
resampling
.
status
.
systematic
if
pruning
...
...
src/conditional_filter_proposal.m
View file @
4921000f
...
...
@@ -51,8 +51,8 @@ ghuu = ReducedForm.ghuu;
ghxu
=
ReducedForm
.
ghxu
;
if
any
(
any
(
isnan
(
ghx
)))
||
any
(
any
(
isnan
(
ghu
)))
||
any
(
any
(
isnan
(
ghxx
)))
||
any
(
any
(
isnan
(
ghuu
)))
||
any
(
any
(
isnan
(
ghxu
)))
||
...
any
(
any
(
isinf
(
ghx
)))
||
any
(
any
(
isinf
(
ghu
)))
||
any
(
any
(
isinf
(
ghxx
)))
||
any
(
any
(
isinf
(
ghuu
)))
||
any
(
any
(
isinf
(
ghxu
)))
...
any
(
any
(
abs
(
ghx
)
>
1e4
))
||
any
(
any
(
abs
(
ghu
)
>
1e4
))
||
any
(
any
(
abs
(
ghxx
)
>
1e4
))
||
any
(
any
(
abs
(
ghuu
)
>
1e4
))
||
any
(
any
(
abs
(
ghxu
)
>
1e4
))
any
(
any
(
isinf
(
ghx
)))
||
any
(
any
(
isinf
(
ghu
)))
||
any
(
any
(
isinf
(
ghxx
)))
||
any
(
any
(
isinf
(
ghuu
)))
||
any
(
any
(
isinf
(
ghxu
)))
...
any
(
any
(
abs
(
ghx
)
>
1e4
))
||
any
(
any
(
abs
(
ghu
)
>
1e4
))
||
any
(
any
(
abs
(
ghxx
)
>
1e4
))
||
any
(
any
(
abs
(
ghuu
)
>
1e4
))
||
any
(
any
(
abs
(
ghxu
)
>
1e4
))
ghx
ghu
ghxx
...
...
@@ -106,7 +106,7 @@ if ParticleOptions.proposal_approximation.cubature || ParticleOptions.proposal_a
Error
=
obs
-
PredictedObservedMean
;
StateVectorMean
=
PredictedStateMean
+
(
CovarianceObservedStateSquareRoot
/
PredictedObservedVarianceSquareRoot
)
*
Error
;
if
ParticleOptions
.
cpf_weights_method
.
amisanotristani
Weights
=
SampleWeights
.*
probability2
(
zeros
(
number_of_observed_variables
,
1
),
PredictedObservedVarianceSquareRoot
,
Error
)
;
Weights
=
SampleWeights
.*
probability2
(
zeros
(
number_of_observed_variables
,
1
),
PredictedObservedVarianceSquareRoot
,
Error
)
;
end
else
dState
=
bsxfun
(
@
minus
,
tmp
(
mf0
,:),
PredictedStateMean
);
...
...
@@ -120,15 +120,15 @@ else
StateVectorVariance
=
PredictedStateVariance
-
KalmanFilterGain
*
PredictedObservedVariance
*
KalmanFilterGain
'
;
StateVectorVarianceSquareRoot
=
chol
(
StateVectorVariance
+
eye
(
number_of_state_variables
)
*
1e-6
)
'
;
if
ParticleOptions
.
cpf_weights_method
.
amisanotristani
Weights
=
SampleWeights
.*
probability2
(
zeros
(
number_of_observed_variables
,
1
),
chol
(
PredictedObservedVariance
)
'
,
Error
)
;
Weights
=
SampleWeights
.*
probability2
(
zeros
(
number_of_observed_variables
,
1
),
chol
(
PredictedObservedVariance
)
'
,
Error
)
;
end
end
PredictedStateVarianceSquareRoot
=
chol
(
PredictedStateVariance
+
eye
(
number_of_state_variables
)
*
1e-6
)
'
;
ProposalStateVector
=
StateVectorVarianceSquareRoot
*
randn
(
size
(
StateVectorVarianceSquareRoot
,
2
),
1
)
+
StateVectorMean
;
if
ParticleOptions
.
cpf_weights_method
.
murrayjonesparslow
Prior
=
probability2
(
PredictedStateMean
,
PredictedStateVarianceSquareRoot
,
ProposalStateVector
)
;
Posterior
=
probability2
(
StateVectorMean
,
StateVectorVarianceSquareRoot
,
ProposalStateVector
)
;
Likelihood
=
probability2
(
obs
,
H_lower_triangular_cholesky
,
measurement_equations
(
ProposalStateVector
,
ReducedForm
,
ThreadsOptions
))
;
Prior
=
probability2
(
PredictedStateMean
,
PredictedStateVarianceSquareRoot
,
ProposalStateVector
)
;
Posterior
=
probability2
(
StateVectorMean
,
StateVectorVarianceSquareRoot
,
ProposalStateVector
)
;
Likelihood
=
probability2
(
obs
,
H_lower_triangular_cholesky
,
measurement_equations
(
ProposalStateVector
,
ReducedForm
,
ThreadsOptions
))
;
Weights
=
SampleWeights
.*
Likelihood
.*
(
Prior
.
/
Posterior
)
;
end
src/conditional_particle_filter.m
View file @
4921000f
function
[
LIK
,
lik
]
=
conditional_particle_filter
(
ReducedForm
,
Y
,
start
,
ParticleOptions
,
ThreadsOptions
)
%
%
% Evaluates the likelihood of a non-linear model with a particle filter
% - the proposal is built using the Kalman updating step for each particle.
% - we need draws in the errors distributions
% Whether we use Monte-Carlo draws from a multivariate gaussian distribution
% as in Amisano & Tristani (JEDC 2010).
% Whether we use multidimensional Gaussian sparse grids approximations:
% - a univariate Kronrod-Paterson Gaussian quadrature combined by the Smolyak
% operator (ref: Winschel & Kratzig, 2010).
% - the proposal is built using the Kalman updating step for each particle.
% - we need draws in the errors distributions
% Whether we use Monte-Carlo draws from a multivariate gaussian distribution
% as in Amisano & Tristani (JEDC 2010).
% Whether we use multidimensional Gaussian sparse grids approximations:
% - a univariate Kronrod-Paterson Gaussian quadrature combined by the Smolyak
% operator (ref: Winschel & Kratzig, 2010).
% - a spherical-radial cubature (ref: Arasaratnam & Haykin, 2009a,2009b).
% - a scaled unscented transform cubature (ref: Julier & Uhlmann 1997, van der
% - a scaled unscented transform cubature (ref: Julier & Uhlmann 1997, van der
% Merwe & Wan 2003).
%
% Pros:
% - Allows using current observable information in the proposal
% - The use of sparse grids Gaussian approximation is much faster than the Monte-Carlo approach
% Cons:
% - The use of the Kalman updating step may biais the proposal distribution since
%
% Pros:
% - Allows using current observable information in the proposal
% - The use of sparse grids Gaussian approximation is much faster than the Monte-Carlo approach
% Cons:
% - The use of the Kalman updating step may biais the proposal distribution since
% it has been derived in a linear context and is implemented in a nonlinear
% context. That is why particle resampling is performed.
% context. That is why particle resampling is performed.
%
% INPUTS
% reduced_form_model [structure] Matlab's structure describing the reduced form model.
...
...
@@ -58,8 +58,8 @@ function [LIK,lik] = conditional_particle_filter(ReducedForm,Y,start,ParticleOpt
% stephane DOT adjemian AT univ DASH lemans DOT fr
persistent
init_flag
mf1
persistent
number_of_particles
persistent
sample_size
number_of_observed_variables
persistent
number_of_particles
persistent
sample_size
number_of_observed_variables
% Set default
if
isempty
(
start
)
...
...
@@ -82,14 +82,14 @@ if isempty(H)
H
=
0
;
H_lower_triangular_cholesky
=
0
;
else
H_lower_triangular_cholesky
=
chol
(
H
)
'
;
H_lower_triangular_cholesky
=
chol
(
H
)
'
;
end
% Get initial condition for the state vector.
StateVectorMean
=
ReducedForm
.
StateVectorMean
;
StateVectorVarianceSquareRoot
=
chol
(
ReducedForm
.
StateVectorVariance
)
'
;
state_variance_rank
=
size
(
StateVectorVarianceSquareRoot
,
2
);
Q_lower_triangular_cholesky
=
chol
(
Q
)
'
;
Q_lower_triangular_cholesky
=
chol
(
Q
)
'
;
% Set seed for randn().
set_dynare_seed
(
'default'
);
...
...
@@ -102,13 +102,13 @@ ks = 0 ;
StateParticles
=
bsxfun
(
@
plus
,
StateVectorVarianceSquareRoot
*
randn
(
state_variance_rank
,
number_of_particles
),
StateVectorMean
);
SampleWeights
=
ones
(
1
,
number_of_particles
)/
number_of_particles
;
for
t
=
1
:
sample_size
for
i
=
1
:
number_of_particles
[
StateParticles
(:,
i
),
SampleWeights
(
i
)]
=
...
conditional_filter_proposal
(
ReducedForm
,
Y
(:,
t
),
StateParticles
(:,
i
),
SampleWeights
(
i
),
Q_lower_triangular_cholesky
,
H_lower_triangular_cholesky
,
H
,
ParticleOptions
,
ThreadsOptions
,
normconst2
)
;
for
i
=
1
:
number_of_particles
[
StateParticles
(:,
i
),
SampleWeights
(
i
)]
=
...
conditional_filter_proposal
(
ReducedForm
,
Y
(:,
t
),
StateParticles
(:,
i
),
SampleWeights
(
i
),
Q_lower_triangular_cholesky
,
H_lower_triangular_cholesky
,
H
,
ParticleOptions
,
ThreadsOptions
,
normconst2
)
;
end
SumSampleWeights
=
sum
(
SampleWeights
)
;
lik
(
t
)
=
log
(
SumSampleWeights
)
;
SampleWeights
=
SampleWeights
.
/
SumSampleWeights
;
lik
(
t
)
=
log
(
SumSampleWeights
)
;
SampleWeights
=
SampleWeights
.
/
SumSampleWeights
;
if
(
ParticleOptions
.
resampling
.
status
.
generic
&&
neff
(
SampleWeights
)
<
ParticleOptions
.
resampling
.
threshold
*
sample_size
)
||
ParticleOptions
.
resampling
.
status
.
systematic
ks
=
ks
+
1
;
StateParticles
=
resample
(
StateParticles
',SampleWeights'
,
ParticleOptions
)
'
;
...
...
src/fit_gaussian_mixture.m
View file @
4921000f
function
[
StateMu
,
StateSqrtP
,
StateWeights
]
=
fit_gaussian_mixture
(
X
,
X_weights
,
StateMu
,
StateSqrtP
,
StateWeights
,
crit
,
niters
,
check
)
function
[
StateMu
,
StateSqrtP
,
StateWeights
]
=
fit_gaussian_mixture
(
X
,
X_weights
,
StateMu
,
StateSqrtP
,
StateWeights
,
crit
,
niters
,
check
)
% Copyright (C) 2013 Dynare Team
%
...
...
@@ -17,36 +17,35 @@ function [StateMu,StateSqrtP,StateWeights] = fit_gaussian_mixture(X,X_weights,St
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
[
dim
,
Ndata
]
=
size
(
X
);
[
dim
,
Ndata
]
=
size
(
X
);
M
=
size
(
StateMu
,
2
)
;
if
check
% Ensure that covariances don't collapse
MIN_COVAR_SQRT
=
sqrt
(
eps
);
init_covars
=
StateSqrtP
;
MIN_COVAR_SQRT
=
sqrt
(
eps
);
init_covars
=
StateSqrtP
;
end
eold
=
-
Inf
;
for
n
=
1
:
niters
% Calculate posteriors based on old parameters
[
prior
,
likelihood
,
marginal
,
posterior
]
=
probability3
(
StateMu
,
StateSqrtP
,
StateWeights
,
X
,
X_weights
);
e
=
sum
(
log
(
marginal
));
if
(
n
>
1
&&
abs
((
e
-
eold
)/
eold
)
<
crit
)
return
;
else
eold
=
e
;
end
new_pr
=
(
sum
(
posterior
,
2
))
'
;
StateWeights
=
new_pr
/
Ndata
;
StateMu
=
bsxfun
(
@
rdivide
,(
posterior
*
X
')'
,
new_pr
);
for
j
=
1
:
M
diffs
=
bsxfun
(
@
minus
,
X
,
StateMu
(:,
j
));
tpost
=
(
1
/
sqrt
(
new_pr
(
j
)))
*
sqrt
(
posterior
(
j
,:));
diffs
=
bsxfun
(
@
times
,
diffs
,
tpost
);
[
foo
,
tcov
]
=
qr2
(
diffs
'
,
0
);
StateSqrtP
(:,:,
j
)
=
tcov
'
;
if
check
if
min
(
abs
(
diag
(
StateSqrtP
(:,:,
j
))))
<
MIN_COVAR_SQRT
StateSqrtP
(:,:,
j
)
=
init_covars
(:,:,
j
);
end
% Calculate posteriors based on old parameters
[
prior
,
likelihood
,
marginal
,
posterior
]
=
probability3
(
StateMu
,
StateSqrtP
,
StateWeights
,
X
,
X_weights
);
e
=
sum
(
log
(
marginal
));
if
(
n
>
1
&&
abs
((
e
-
eold
)/
eold
)
<
crit
)
return
;
else
eold
=
e
;
end
end
end
new_pr
=
(
sum
(
posterior
,
2
))
'
;
StateWeights
=
new_pr
/
Ndata
;
StateMu
=
bsxfun
(
@
rdivide
,(
posterior
*
X
')'
,
new_pr
);
for
j
=
1
:
M
diffs
=
bsxfun
(
@
minus
,
X
,
StateMu
(:,
j
));
tpost
=
(
1
/
sqrt
(
new_pr
(
j
)))
*
sqrt
(
posterior
(
j
,:));
diffs
=
bsxfun
(
@
times
,
diffs
,
tpost
);
[
foo
,
tcov
]
=
qr2
(
diffs
'
,
0
);
StateSqrtP
(:,:,
j
)
=
tcov
'
;
if
check
if
min
(
abs
(
diag
(
StateSqrtP
(:,:,
j
))))
<
MIN_COVAR_SQRT
StateSqrtP
(:,:,
j
)
=
init_covars
(:,:,
j
);
end
end
end
end
src/gaussian_densities.m
View file @
4921000f
function
IncrementalWeights
=
gaussian_densities
(
obs
,
mut_t
,
sqr_Pss_t_t
,
st_t_1
,
sqr_Pss_t_t_1
,
particles
,
H
,
normconst
,
weigths1
,
weigths2
,
ReducedForm
,
ThreadsOptions
)
%
% Elements to calculate the importance sampling ratio
% Elements to calculate the importance sampling ratio
%
% INPUTS
% reduced_form_model [structure] Matlab's structure describing the reduced form model.
...
...
@@ -36,11 +36,11 @@ function IncrementalWeights = gaussian_densities(obs,mut_t,sqr_Pss_t_t,st_t_1,sq
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
% proposal density
proposal
=
probability2
(
mut_t
,
sqr_Pss_t_t
,
particles
)
;
% prior density
prior
=
probability2
(
st_t_1
,
sqr_Pss_t_t_1
,
particles
)
;
% likelihood
% proposal density
proposal
=
probability2
(
mut_t
,
sqr_Pss_t_t
,
particles
)
;
% prior density
prior
=
probability2
(
st_t_1
,
sqr_Pss_t_t_1
,
particles
)
;
% likelihood
yt_t_1_i
=
measurement_equations
(
particles
,
ReducedForm
,
ThreadsOptions
)
;
eta_t_i
=
bsxfun
(
@
minus
,
obs
,
yt_t_1_i
)
'
;
yt_t_1
=
sum
(
yt_t_1_i
*
weigths1
,
2
)
;
...
...
@@ -48,5 +48,5 @@ tmp = bsxfun(@minus,yt_t_1_i,yt_t_1) ;
Pyy
=
bsxfun
(
@
times
,
weigths2
',tmp)*tmp'
+
H
;
sqr_det
=
sqrt
(
det
(
Pyy
))
;
foo
=
(
eta_t_i
/
Pyy
)
.*
eta_t_i
;
likelihood
=
exp
(
-
0.5
*
sum
(
foo
,
2
))/(
normconst
*
sqr_det
)
+
1e-99
;
likelihood
=
exp
(
-
0.5
*
sum
(
foo
,
2
))/(
normconst
*
sqr_det
)
+
1e-99
;
IncrementalWeights
=
likelihood
.*
prior
.
/
proposal
;
src/gaussian_filter.m
View file @
4921000f
...
...
@@ -112,18 +112,18 @@ for t=1:sample_size
if
ParticleOptions
.
distribution_approximation
.
cubature
||
ParticleOptions
.
distribution_approximation
.
unscented
StateParticles
=
bsxfun
(
@
plus
,
StateVectorMean
,
StateVectorVarianceSquareRoot
*
nodes2
'
)
;
IncrementalWeights
=
...
gaussian_densities
(
Y
(:,
t
),
StateVectorMean
,
...
StateVectorVarianceSquareRoot
,
PredictedStateMean
,
...
PredictedStateVarianceSquareRoot
,
StateParticles
,
H
,
const_lik
,
...
weights2
,
weights_c2
,
ReducedForm
,
ThreadsOptions
)
;
gaussian_densities
(
Y
(:,
t
),
StateVectorMean
,
...
StateVectorVarianceSquareRoot
,
PredictedStateMean
,
...
PredictedStateVarianceSquareRoot
,
StateParticles
,
H
,
const_lik
,
...
weights2
,
weights_c2
,
ReducedForm
,
ThreadsOptions
)
;
SampleWeights
=
weights2
.*
IncrementalWeights
;
else
else
StateParticles
=
bsxfun
(
@
plus
,
StateVectorVarianceSquareRoot
*
randn
(
state_variance_rank
,
number_of_particles
),
StateVectorMean
)
;
IncrementalWeights
=
...
gaussian_densities
(
Y
(:,
t
),
StateVectorMean
,
...
StateVectorVarianceSquareRoot
,
PredictedStateMean
,
...
PredictedStateVarianceSquareRoot
,
StateParticles
,
H
,
const_lik
,
...
1
/
number_of_particles
,
1
/
number_of_particles
,
ReducedForm
,
ThreadsOptions
)
;
gaussian_densities
(
Y
(:,
t
),
StateVectorMean
,
...
StateVectorVarianceSquareRoot
,
PredictedStateMean
,
...
PredictedStateVarianceSquareRoot
,
StateParticles
,
H
,
const_lik
,
...
1
/
number_of_particles
,
1
/
number_of_particles
,
ReducedForm
,
ThreadsOptions
)
;
SampleWeights
=
IncrementalWeights
/
number_of_particles
;
end
SampleWeights
=
SampleWeights
+
1e-6
*
ones
(
size
(
SampleWeights
,
1
),
1
)
;
...
...
@@ -132,9 +132,9 @@ for t=1:sample_size
SampleWeights
=
SampleWeights
.
/
SumSampleWeights
;
if
not
(
ParticleOptions
.
distribution_approximation
.
cubature
||
ParticleOptions
.
distribution_approximation
.
unscented
)
if
(
ParticleOptions
.
resampling
.
status
.
generic
&&
neff
(
SampleWeights
)
<
ParticleOptions
.
resampling
.
threshold
*
sample_size
)
||
ParticleOptions
.
resampling
.
status
.
systematic
StateParticles
=
resample
(
StateParticles
',SampleWeights,ParticleOptions)'
;
StateParticles
=
resample
(
StateParticles
',SampleWeights,ParticleOptions)'
;
SampleWeights
=
ones
(
number_of_particles
,
1
)/
number_of_particles
;
end
end
end
StateVectorMean
=
StateParticles
*
SampleWeights
;
temp
=
bsxfun
(
@
minus
,
StateParticles
,
StateVectorMean
)
;
...
...
src/gaussian_filter_bank.m
View file @
4921000f
...
...
@@ -49,8 +49,8 @@ ghuu = ReducedForm.ghuu;
ghxu
=
ReducedForm
.
ghxu
;
if
any
(
any
(
isnan
(
ghx
)))
||
any
(
any
(
isnan
(
ghu
)))
||
any
(
any
(
isnan
(
ghxx
)))
||
any
(
any
(
isnan
(
ghuu
)))
||
any
(
any
(
isnan
(
ghxu
)))
||
...
any
(
any
(
isinf
(
ghx
)))
||
any
(
any
(
isinf
(
ghu
)))
||
any
(
any
(
isinf
(
ghxx
)))
||
any
(
any
(
isinf
(
ghuu
)))
||
any
(
any
(
isinf
(
ghxu
)))
...
any
(
any
(
abs
(
ghx
)
>
1e4
))
||
any
(
any
(
abs
(
ghu
)
>
1e4
))
||
any
(
any
(
abs
(
ghxx
)
>
1e4
))
||
any
(
any
(
abs
(
ghuu
)
>
1e4
))
||
any
(
any
(
abs
(
ghxu
)
>
1e4
))
any
(
any
(
isinf
(
ghx
)))
||
any
(
any
(
isinf
(
ghu
)))
||
any
(
any
(
isinf
(
ghxx
)))
||
any
(
any
(
isinf
(
ghuu
)))
||
any
(
any
(
isinf
(
ghxu
)))
...
any
(
any
(
abs
(
ghx
)
>
1e4
))
||
any
(
any
(
abs
(
ghu
)
>
1e4
))
||
any
(
any
(
abs
(
ghxx
)
>
1e4
))
||
any
(
any
(
abs
(
ghuu
)
>
1e4
))
||
any
(
any
(
abs
(
ghxu
)
>
1e4
))
ghx
ghu
ghxx
...
...
src/gaussian_mixture_densities.m
View file @
4921000f
function
IncrementalWeights
=
gaussian_mixture_densities
(
obs
,
StateMuPrior
,
StateSqrtPPrior
,
StateWeightsPrior
,
...
StateMuPost
,
StateSqrtPPost
,
StateWeightsPost
,
...
StateParticles
,
H
,
normconst
,
weigths1
,
weigths2
,
ReducedForm
,
ThreadsOptions
)
StateMuPost
,
StateSqrtPPost
,
StateWeightsPost
,
...
StateParticles
,
H
,
normconst
,
weigths1
,
weigths2
,
ReducedForm
,
ThreadsOptions
)
%
% Elements to calculate the importance sampling ratio
% Elements to calculate the importance sampling ratio
%
% INPUTS
% reduced_form_model [structure] Matlab's structure describing the reduced form model.
...
...
@@ -38,11 +38,11 @@ function IncrementalWeights = gaussian_mixture_densities(obs,StateMuPrior,State
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
% Compute the density of particles under the prior distribution
[
ras
,
ras
,
prior
]
=
probability
(
StateMuPrior
,
StateSqrtPPrior
,
StateWeightsPrior
,
StateParticles
)
;
% Compute the density of particles under the prior distribution
[
ras
,
ras
,
prior
]
=
probability
(
StateMuPrior
,
StateSqrtPPrior
,
StateWeightsPrior
,
StateParticles
)
;
prior
=
prior
'
;
% Compute the density of particles under the proposal distribution
[
ras
,
ras
,
proposal
]
=
probability
(
StateMuPost
,
StateSqrtPPost
,
StateWeightsPost
,
StateParticles
)
;
% Compute the density of particles under the proposal distribution
[
ras
,
ras
,
proposal
]
=
probability
(
StateMuPost
,
StateSqrtPPost
,
StateWeightsPost
,
StateParticles
)
;
proposal
=
proposal
'
;
% Compute the density of the current observation conditionally to each particle
yt_t_1_i
=
measurement_equations
(
StateParticles
,
ReducedForm
,
ThreadsOptions
)
;
...
...
@@ -52,6 +52,5 @@ tmp = bsxfun(@minus,yt_t_1_i,yt_t_1) ;
Pyy
=
bsxfun
(
@
times
,
weigths2
',tmp)*tmp'
+
H
;
sqr_det
=
sqrt
(
det
(
Pyy
))
;
foo
=
(
eta_t_i
/
Pyy
)
.*
eta_t_i
;
likelihood
=
exp
(
-
0.5
*
sum
(
foo
,
2
))/(
normconst
*
sqr_det
)
+
1e-99
;
likelihood
=
exp
(
-
0.5
*
sum
(
foo
,
2
))/(
normconst
*
sqr_det
)
+
1e-99
;
IncrementalWeights
=
likelihood
.*
prior
.
/
proposal
;
src/gaussian_mixture_filter.m
View file @
4921000f
...
...
@@ -63,15 +63,15 @@ end
% Set persistent variables.
if
isempty
(
init_flag
)
mf0
=
ReducedForm
.
mf0
;
mf1
=
ReducedForm
.
mf1
;
sample_size
=
size
(
Y
,
2
);
number_of_state_variables
=
length
(
mf0
);
number_of_observed_variables
=
length
(
mf1
);
number_of_structural_innovations
=
length
(
ReducedForm
.
Q
);
G
=
ParticleOptions
.
mixture_state_variables
;
% number of GM components in state
number_of_particles
=
ParticleOptions
.
number_of_particles
;
init_flag
=
1
;
mf0
=
ReducedForm
.
mf0
;
mf1
=
ReducedForm
.
mf1
;
sample_size
=
size
(
Y
,
2
);
number_of_state_variables
=
length
(
mf0
);
number_of_observed_variables
=
length
(
mf1
);
number_of_structural_innovations
=
length
(
ReducedForm
.
Q
);
G
=
ParticleOptions
.
mixture_state_variables
;
% number of GM components in state
number_of_particles
=
ParticleOptions
.
number_of_particles
;
init_flag
=
1
;
end
% compute gaussian quadrature nodes and weights on states and shocks
...
...
@@ -104,14 +104,14 @@ else
end
Q_lower_triangular_cholesky
=
reduced_rank_cholesky
(
Q
)
'
;
% Initialize mixtures
% Initialize mixtures
StateWeights
=
ones
(
1
,
G
)/
G
;
StateMu
=
ReducedForm
.
StateVectorMean
;
StateSqrtP
=
zeros
(
number_of_state_variables
,
number_of_state_variables
,
G
)
;
temp
=
reduced_rank_cholesky
(
ReducedForm
.
StateVectorVariance
)
'
;
StateMu
=
bsxfun
(
@
plus
,
StateMu
,
bsxfun
(
@
times
,
diag
(
temp
),(
-
(
G
-
1
)/
2
:
1
:(
G
-
1
)/
2
))/
10
)
;
for
g
=
1
:
G
StateSqrtP
(:,:,
g
)
=
temp
/
sqrt
(
G
)
;
StateSqrtP
(:,:,
g
)
=
temp
/
sqrt
(
G
)
;
end
% if ParticleOptions.mixture_structural_shocks==1
...
...
@@ -135,11 +135,11 @@ end
% for i=1:I
% StructuralShocksSqrtP(:,:,i) = Q_lower_triangular_cholesky/sqrt(StructuralShocksWeights(i)) ;
% end
%
% if ParticleOptions.mixture_measurement_shocks==1
%
% if ParticleOptions.mixture_measurement_shocks==1
% ObservationShocksMu = zeros(1,number_of_observed_variables) ;
% ObservationShocksWeights = 1 ;
% else
% else
% if ParticleOptions.proposal_approximation.cubature
% [ObservationShocksMu,ObservationShocksWeights] = spherical_radial_sigma_points(number_of_observed_variables);
% ObservationShocksWeights = ones(size(ObservationShocksMu,1),1)*ObservationShocksWeights;
...
...
@@ -150,7 +150,7 @@ end
% error('Estimation: This approximation for the proposal is not implemented or unknown!')
% end
% end
% end
% end
% J = size(ObservationShocksWeights,1) ;
% ObservationShocksMu = H_lower_triangular_cholesky*(ObservationShocksMu') ;
% ObservationShocksSqrtP = zeros(number_of_observed_variables,number_of_observed_variables,J) ;
...
...
@@ -180,7 +180,7 @@ elseif ParticleOptions.mixture_structural_shocks==1
StructuralShocksMu
=
Q_lower_triangular_cholesky
*
(
StructuralShocksMu
'
)
;
StructuralShocksSqrtP
=
zeros
(
number_of_structural_innovations
,
number_of_structural_innovations
,
I
)
;
for
i
=
1
:
I
StructuralShocksSqrtP
(:,:,
i
)
=
Q_lower_triangular_cholesky
;
StructuralShocksSqrtP
(:,:,
i
)
=
Q_lower_triangular_cholesky
;
end
else
if
ParticleOptions
.
proposal_approximation
.
cubature
...
...
@@ -197,7 +197,7 @@ else
StructuralShocksMu
=
Q_lower_triangular_cholesky
*
(
StructuralShocksMu
'
)
;
StructuralShocksSqrtP
=
zeros
(
number_of_structural_innovations
,
number_of_structural_innovations
,
I
)
;
for
i
=
1
:
I
StructuralShocksSqrtP
(:,:,
i
)
=
Q_lower_triangular_cholesky
/
sqrt
(
StructuralShocksWeights
(
i
))
;
StructuralShocksSqrtP
(:,:,
i
)
=
Q_lower_triangular_cholesky
/
sqrt
(
StructuralShocksWeights
(
i
))
;
end
end
...
...
@@ -208,14 +208,14 @@ ObservationShocksMu = H_lower_triangular_cholesky*(ObservationShocksMu') ;
ObservationShocksSqrtP
=
zeros
(
number_of_observed_variables
,
number_of_observed_variables
,
J
)
;
ObservationShocksSqrtP
(:,:,
1
)
=
H_lower_triangular_cholesky
;
% if ParticleOptions.mixture_measurement_shocks==0
% if ParticleOptions.mixture_measurement_shocks==0
% ObservationShocksMu = zeros(1,number_of_observed_variables) ;
% ObservationShocksWeights = 1 ;
% J = 1 ;
% ObservationShocksMu = H_lower_triangular_cholesky*(ObservationShocksMu') ;
% ObservationShocksSqrtP = zeros(number_of_observed_variables,number_of_observed_variables,J) ;
% ObservationShocksSqrtP(:,:,1) = H_lower_triangular_cholesky ;
% elseif ParticleOptions.mixture_measurement_shocks==1
% elseif ParticleOptions.mixture_measurement_shocks==1
% if ParticleOptions.proposal_approximation.cubature
% [ObservationShocksMu,ObservationShocksWeights] = spherical_radial_sigma_points(number_of_observed_variables);
% ObservationShocksWeights = ones(size(ObservationShocksMu,1),1)*ObservationShocksWeights;
...
...
@@ -232,7 +232,7 @@ ObservationShocksSqrtP(:,:,1) = H_lower_triangular_cholesky ;
% for j=1:J
% ObservationShocksSqrtP(:,:,j) = H_lower_triangular_cholesky ;
% end
% else
% else
% if ParticleOptions.proposal_approximation.cubature
% [ObservationShocksMu,ObservationShocksWeights] = spherical_radial_sigma_points(number_of_observed_variables);
% ObservationShocksWeights = ones(size(ObservationShocksMu,1),1)*ObservationShocksWeights;
...
...
@@ -277,10 +277,10 @@ for t=1:sample_size
gsecond
=
gprime
+
(
j
-
1
)
*
Gprime
;
[
StateMuPrior
(:,
gprime
),
StateSqrtPPrior
(:,:,
gprime
),
StateWeightsPrior
(
1
,
gprime
),
...
StateMuPost
(:,
gsecond
),
StateSqrtPPost
(:,:,
gsecond
),
StateWeightsPost
(
1
,
gsecond
)]
=
...
gaussian_mixture_filter_bank
(
ReducedForm
,
Y
(:,
t
),
StateMu
(:,
g
),
StateSqrtP
(:,:,
g
),
StateWeights
(
g
),
...
StructuralShocksMu
(:,
i
),
StructuralShocksSqrtP
(:,:,
i
),
StructuralShocksWeights
(
i
),
...
ObservationShocksMu
(:,
j
),
ObservationShocksSqrtP
(:,:,
j
),
ObservationShocksWeights
(
j
),
...
H
,
H_lower_triangular_cholesky
,
const_lik
,
ParticleOptions
,
ThreadsOptions
)
;
gaussian_mixture_filter_bank
(
ReducedForm
,
Y
(:,
t
),
StateMu
(:,
g
),
StateSqrtP
(:,:,
g
),
StateWeights
(
g
),
...
StructuralShocksMu
(:,
i
),
StructuralShocksSqrtP
(:,:,
i
),
StructuralShocksWeights
(
i
),
...
ObservationShocksMu
(:,
j
),
ObservationShocksSqrtP
(:,:,
j
),
ObservationShocksWeights
(
j
),
...
H
,
H_lower_triangular_cholesky
,
const_lik
,
ParticleOptions
,
ThreadsOptions
)
;
end
end
end
...
...
@@ -293,8 +293,8 @@ for t=1:sample_size
for
i
=
1
:
Gsecond
StateParticles
=
bsxfun
(
@
plus
,
StateMuPost
(:,
i
),
StateSqrtPPost
(:,:,
i
)
*
nodes
'
)
;
IncrementalWeights
=
gaussian_mixture_densities
(
Y
(:,
t
),
StateMuPrior
,
StateSqrtPPrior
,
StateWeightsPrior
,
...
StateMuPost
,
StateSqrtPPost
,
StateWeightsPost
,
...
StateParticles
,
H
,
const_lik
,
weights
,
weights_c
,
ReducedForm
,
ThreadsOptions
)
;
StateMuPost
,
StateSqrtPPost
,
StateWeightsPost
,
...
StateParticles
,
H
,
const_lik
,
weights
,
weights_c
,
ReducedForm
,
ThreadsOptions
)
;
SampleWeights
(
i
)
=
sum
(
StateWeightsPost
(
i
)
*
weights
.*
IncrementalWeights
)
;
end
SumSampleWeights
=
sum
(
SampleWeights
)
;
...
...
@@ -311,9 +311,9 @@ for t=1:sample_size
% Sample particle in the proposal distribution, ie the posterior state GM
StateParticles
=
importance_sampling
(
StateMuPost
,
StateSqrtPPost
,
StateWeightsPost
'
,
number_of_particles
)
;
IncrementalWeights
=
gaussian_mixture_densities
(
Y
(:,
t
),
StateMuPrior
,
StateSqrtPPrior
,
StateWeightsPrior
,
...
StateMuPost
,
StateSqrtPPost
,
StateWeightsPost
,
...
StateParticles
,
H
,
const_lik
,
1
/
number_of_particles
,
...
1
/
number_of_particles
,
ReducedForm
,
ThreadsOptions
)
;
StateMuPost
,
StateSqrtPPost
,
StateWeightsPost
,
...
StateParticles
,
H
,
const_lik
,
1
/
number_of_particles
,
...
1
/
number_of_particles
,
ReducedForm
,
ThreadsOptions
)
;
SampleWeights
=
IncrementalWeights
/
number_of_particles
;
SumSampleWeights
=
sum
(
SampleWeights
,
1
)
;
SampleWeights
=
SampleWeights
.
/
SumSampleWeights
;
...
...
src/gaussian_mixture_filter_bank.m
View file @
4921000f
function
[
StateMuPrior
,
StateSqrtPPrior
,
StateWeightsPrior
,
StateMuPost
,
StateSqrtPPost
,
StateWeightsPost
]
=
...
gaussian_mixture_filter_bank
(
ReducedForm
,
obs
,
StateMu
,
StateSqrtP
,
StateWeights
,
...
StructuralShocksMu
,
StructuralShocksSqrtP
,
StructuralShocksWeights
,
...
ObservationShocksMu
,
ObservationShocksSqrtP
,
ObservationShocksWeights
,
...
H
,
H_lower_triangular_cholesky
,
normfactO
,
ParticleOptions
,
ThreadsOptions
)
gaussian_mixture_filter_bank
(
ReducedForm
,
obs
,
StateMu
,
StateSqrtP
,
StateWeights
,
...
StructuralShocksMu
,
StructuralShocksSqrtP
,
StructuralShocksWeights
,
...
ObservationShocksMu
,
ObservationShocksSqrtP
,
ObservationShocksWeights
,
...
H
,
H_lower_triangular_cholesky
,
normfactO
,
ParticleOptions
,
ThreadsOptions
)
%
% Computes the proposal with a gaussian approximation for importance
% sampling
% This proposal is a gaussian distribution calculated à la Kalman
% sampling
% This proposal is a gaussian distribution calculated à la Kalman
%
% INPUTS
% reduced_form_model [structure] Matlab's structure describing the reduced form model.
...
...
@@ -43,8 +43,8 @@ function [StateMuPrior,StateSqrtPPrior,StateWeightsPrior,StateMuPost,StateSqrtPP
persistent
init_flag2
mf0
mf1
%nodes3 weights3 weights_c3
persistent
number_of_state_variables
number_of_observed_variables
persistent
number_of_structural_innovations
persistent
number_of_state_variables
number_of_observed_variables
persistent
number_of_structural_innovations
% Set local state space model (first-order approximation).
ghx
=
ReducedForm
.
ghx
;
...
...
@@ -55,8 +55,8 @@ ghuu = ReducedForm.ghuu;
ghxu
=
ReducedForm
.
ghxu
;
if
any
(
any
(
isnan
(
ghx
)))
||
any
(
any
(
isnan
(
ghu
)))
||
any
(
any
(
isnan
(
ghxx
)))
||
any
(
any
(
isnan
(
ghuu
)))
||
any
(
any
(
isnan
(
ghxu
)))
||
...
any
(
any
(
isinf
(
ghx
)))
||
any
(
any
(
isinf
(
ghu
)))
||
any
(
any
(
isinf
(
ghxx
)))
||
any
(
any
(
isinf
(
ghuu
)))
||
any
(
any
(
isinf
(
ghxu
)))
...
any
(
any
(
abs
(
ghx
)
>
1e4
))
||
any
(
any
(
abs
(
ghu
)
>
1e4
))
||
any
(
any
(
abs
(
ghxx
)
>
1e4
))
||
any
(
any
(
abs
(
ghuu
)
>
1e4
))
||
any
(
any
(
abs
(
ghxu
)
>
1e4
))
any
(
any
(
isinf
(
ghx
)))
||
any
(
any
(
isinf
(
ghu
)))
||
any
(
any
(
isinf
(
ghxx
)))
||
any
(
any
(
isinf
(
ghuu
)))
||
any
(
any
(
isinf
(
ghxu
)))
...
any
(
any
(
abs
(
ghx
)
>
1e4
))
||
any
(
any
(
abs
(
ghu
)
>
1e4
))
||
any
(
any
(
abs
(
ghxx
)
>
1e4
))
||
any
(
any
(
abs
(
ghuu
)
>
1e4
))
||
any
(
any
(
abs
(
ghxu
)
>
1e4
))
ghx
ghu
ghxx
...
...
src/importance_sampling.m
View file @
4921000f
...
...
@@ -17,13 +17,12 @@ function State_Particles = importance_sampling(StateMuPost,StateSqrtPPost,StateW
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
[
Xdim
,
Gsecond
]
=
size
(
StateMuPost
)
;
[
Xdim
,
Gsecond
]
=
size
(
StateMuPost
)
;
u
=
rand
(
numP
,
1
);
[
Nc
,
comp
]
=
histc
(
u
,
cumsum
([
0
;
StateWeightsPost
]));
[
Nc
,
comp
]
=
histc
(
u
,
cumsum
([
0
;
StateWeightsPost
]));
State_Particles
=
zeros
(
Xdim
,
numP
);
for
k
=
1
:
Gsecond
idx
=
bsxfun
(
@
eq
,
comp
,
k
*
ones
(
size
(
comp
)))
;
State_Particles
(:,
idx
)
=
StateSqrtPPost
(:,:,
k
)
*
randn
(
Xdim
,
Nc
(
k
));
idx
=
bsxfun
(
@
eq
,
comp
,
k
*
ones
(
size
(
comp
)))
;
State_Particles
(:,
idx
)
=
StateSqrtPPost
(:,:,
k
)
*
randn
(
Xdim
,
Nc
(
k
));
end
State_Particles
=
State_Particles
+
StateMuPost
(:,
comp
);
State_Particles
=
State_Particles
+
StateMuPost
(:,
comp
);
src/measurement_equations.m
View file @
4921000f
function
measure
=
measurement_equations
(
StateVectors
,
ReducedForm
,
ThreadsOptions
)
function
measure
=
measurement_equations
(
StateVectors
,
ReducedForm
,
ThreadsOptions
)
% Copyright (C) 2013 Dynare Team
%
...
...
@@ -28,6 +28,3 @@ state_variables_steady_state = ReducedForm.state_variables_steady_state;
number_of_structural_innovations
=
length
(
ReducedForm
.
Q
);
yhat
=
bsxfun
(
@
minus
,
StateVectors
,
state_variables_steady_state
)
;
measure
=
local_state_space_iteration_2
(
yhat
,
zeros
(
number_of_structural_innovations
,
size
(
yhat
,
2
)),
ghx
,
ghu
,
constant
,
ghxx
,
ghuu
,
ghxu
,
ThreadsOptions
.
local_state_space_iteration_2
);