gaussian_mixture_filter_bank.m 7.03 KB
Newer Older
 Stéphane Adjemian committed Nov 10, 2014 1 ``````function [StateMuPrior,StateSqrtPPrior,StateWeightsPrior,StateMuPost,StateSqrtPPost,StateWeightsPost] =... `````` Stéphane Adjemian committed May 18, 2017 2 3 4 5 `````` gaussian_mixture_filter_bank(ReducedForm,obs,StateMu,StateSqrtP,StateWeights,... StructuralShocksMu,StructuralShocksSqrtP,StructuralShocksWeights,... ObservationShocksMu,ObservationShocksSqrtP,ObservationShocksWeights,... H,H_lower_triangular_cholesky,normfactO,ParticleOptions,ThreadsOptions) `````` Stéphane Adjemian committed Nov 10, 2014 6 7 ``````% % Computes the proposal with a gaussian approximation for importance `````` Stéphane Adjemian committed May 18, 2017 8 9 ``````% sampling % This proposal is a gaussian distribution calculated à la Kalman `````` Stéphane Adjemian committed Nov 10, 2014 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 ``````% % INPUTS % reduced_form_model [structure] Matlab's structure describing the reduced form model. % reduced_form_model.measurement.H [double] (pp x pp) variance matrix of measurement errors. % reduced_form_model.state.Q [double] (qq x qq) variance matrix of state errors. % reduced_form_model.state.dr [structure] output of resol.m. % Y [double] pp*smpl matrix of (detrended) data, where pp is the maximum number of observed variables. % % OUTPUTS % LIK [double] scalar, likelihood % lik [double] vector, density of observations in each period. % % REFERENCES % % NOTES % The vector "lik" is used to evaluate the jacobian of the likelihood. % Copyright (C) 2009-2012 Dynare Team % % This file is part of Dynare. % % Dynare is free software: you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by % the Free Software Foundation, either version 3 of the License, or % (at your option) any later version. % % Dynare is distributed in the hope that it will be useful, % but WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details. % % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . persistent init_flag2 mf0 mf1 %nodes3 weights3 weights_c3 `````` Stéphane Adjemian committed May 18, 2017 46 47 ``````persistent number_of_state_variables number_of_observed_variables persistent number_of_structural_innovations `````` Stéphane Adjemian committed Nov 10, 2014 48 49 50 51 52 53 54 55 56 57 `````` % Set local state space model (first-order approximation). ghx = ReducedForm.ghx; ghu = ReducedForm.ghu; % Set local state space model (second-order approximation). ghxx = ReducedForm.ghxx; 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))) || ... `````` Stéphane Adjemian committed May 18, 2017 58 59 `````` 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)) `````` Stéphane Adjemian committed Nov 10, 2014 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 `````` ghx ghu ghxx ghuu ghxu end constant = ReducedForm.constant; state_variables_steady_state = ReducedForm.state_variables_steady_state; % Set persistent variables. if isempty(init_flag2) mf0 = ReducedForm.mf0; mf1 = ReducedForm.mf1; number_of_state_variables = length(mf0); number_of_observed_variables = length(mf1); number_of_structural_innovations = length(ReducedForm.Q); init_flag2 = 1; end numb = number_of_state_variables+number_of_structural_innovations ; `````` Stéphane Adjemian committed Dec 15, 2014 82 ``````if ParticleOptions.proposal_approximation.cubature `````` Stéphane Adjemian committed Nov 10, 2014 83 84 `````` [nodes3,weights3] = spherical_radial_sigma_points(numb); weights_c3 = weights3; `````` Stéphane Adjemian committed Dec 15, 2014 85 86 ``````elseif ParticleOptions.proposal_approximation.unscented [nodes3,weights3,weights_c3] = unscented_sigma_points(numb,ParticleOptions); `````` Stéphane Adjemian committed Nov 10, 2014 87 88 89 90 91 92 93 ``````else error('Estimation: This approximation for the proposal is not implemented or unknown!') end epsilon = bsxfun(@plus,StructuralShocksSqrtP*nodes3(:,number_of_state_variables+1:number_of_state_variables+number_of_structural_innovations)',StructuralShocksMu) ; StateVectors = bsxfun(@plus,StateSqrtP*nodes3(:,1:number_of_state_variables)',StateMu); yhat = bsxfun(@minus,StateVectors,state_variables_steady_state); `````` Stéphane Adjemian committed Dec 15, 2014 94 ``````tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2); `````` Stéphane Adjemian committed Nov 10, 2014 95 96 97 ``````PredictedStateMean = tmp(mf0,:)*weights3; PredictedObservedMean = tmp(mf1,:)*weights3; `````` Stéphane Adjemian committed Dec 15, 2014 98 ``````if ParticleOptions.proposal_approximation.cubature `````` Stéphane Adjemian committed Nov 10, 2014 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 `````` PredictedStateMean = sum(PredictedStateMean,2); PredictedObservedMean = sum(PredictedObservedMean,2); dState = (bsxfun(@minus,tmp(mf0,:),PredictedStateMean)').*sqrt(weights3); dObserved = (bsxfun(@minus,tmp(mf1,:),PredictedObservedMean)').*sqrt(weights3); PredictedStateVariance = dState'*dState; big_mat = [dObserved dState ; [H_lower_triangular_cholesky zeros(number_of_observed_variables,number_of_state_variables)] ]; [mat1,mat] = qr2(big_mat,0); mat = mat'; clear('mat1'); PredictedObservedVarianceSquareRoot = mat(1:number_of_observed_variables,1:number_of_observed_variables); CovarianceObservedStateSquareRoot = mat(number_of_observed_variables+(1:number_of_state_variables),1:number_of_observed_variables); StateVectorVarianceSquareRoot = mat(number_of_observed_variables+(1:number_of_state_variables),number_of_observed_variables+(1:number_of_state_variables)); iPredictedObservedVarianceSquareRoot = inv(PredictedObservedVarianceSquareRoot); iPredictedObservedVariance = iPredictedObservedVarianceSquareRoot'*iPredictedObservedVarianceSquareRoot; sqrdet = 1/sqrt(det(iPredictedObservedVariance)); PredictionError = obs - PredictedObservedMean; StateVectorMean = PredictedStateMean + CovarianceObservedStateSquareRoot*iPredictedObservedVarianceSquareRoot*PredictionError; else dState = bsxfun(@minus,tmp(mf0,:),PredictedStateMean); dObserved = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean); PredictedStateVariance = dState*diag(weights_c3)*dState'; PredictedObservedVariance = dObserved*diag(weights_c3)*dObserved' + H; PredictedStateAndObservedCovariance = dState*diag(weights_c3)*dObserved'; sqrdet = sqrt(det(PredictedObservedVariance)) ; iPredictedObservedVariance = inv(PredictedObservedVariance); PredictionError = obs - PredictedObservedMean; KalmanFilterGain = PredictedStateAndObservedCovariance*iPredictedObservedVariance; StateVectorMean = PredictedStateMean + KalmanFilterGain*PredictionError; StateVectorVariance = PredictedStateVariance - KalmanFilterGain*PredictedObservedVariance*KalmanFilterGain'; StateVectorVariance = .5*(StateVectorVariance+StateVectorVariance'); StateVectorVarianceSquareRoot = reduced_rank_cholesky(StateVectorVariance)'; end data_lik_GM_g = exp(-0.5*PredictionError'*iPredictedObservedVariance*PredictionError)/abs(normfactO*sqrdet) + 1e-99; StateMuPrior = PredictedStateMean ; StateSqrtPPrior = reduced_rank_cholesky(PredictedStateVariance)'; StateWeightsPrior = StateWeights*StructuralShocksWeights; StateMuPost = StateVectorMean; StateSqrtPPost = StateVectorVarianceSquareRoot; StateWeightsPost = StateWeightsPrior*ObservationShocksWeights*data_lik_GM_g ;``````