stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge remote-tracking branch 'ratto/master'

time-shift
Sébastien Villemot 2012-06-08 18:24:18 +02:00
parent 1f9cea669a bf88b6e93d
commit 129553579a
21 changed files with 243 additions and 166 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
matlab/AHessian.m
View File

@ -50,7 +50,7 @@ end
% DOm = DR(:,:,ii)*Q*transpose(R) + R*DQ(:,:,ii)*transpose(R) + R*Q*transpose(DR(:,:,ii));
% end
while notsteady & t<smpl
while notsteady && t<smpl
t = t+1;
v = Y(:,t)-a(mf);
F = P(mf,mf) + H;
@ -112,7 +112,7 @@ end
a = T*(a+K*v);
lik(t) = transpose(v)*iF*v;
end
AHess = AHess + .5*(smpl+t0-1)*(vecDPmf' * kron(iF,iF) * vecDPmf);
AHess = AHess + .5*(smpl-t0+1)*(vecDPmf' * kron(iF,iF) * vecDPmf);
if nargout > 1
for ii = 1:k
% DLIK(ii,1) = DLIK(ii,1) + (smpl-t0+1)*trace( iF*DF(:,:,ii) );

32
matlab/dsge_likelihood.m
View File

@ -332,7 +332,7 @@ if (kalman_algo == 2) || (kalman_algo == 4)
else
if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
H = diag(H);
mmm = mm;
mmm = mm;
else
Z = [Z, eye(pp)];
T = blkdiag(T,zeros(pp));
@ -381,7 +381,7 @@ switch DynareOptions.lik_init
% Use standard kalman filter except if the univariate filter is explicitely choosen.
if kalman_algo == 0
kalman_algo = 3;
elseif ~((kalman_algo == 3) || (kalman_algo == 4))
elseif ~((kalman_algo == 3) || (kalman_algo == 4))
error(['diffuse filter: options_.kalman_algo can only be equal ' ...
'to 0 (default), 3 or 4'])
end
@ -455,7 +455,7 @@ switch DynareOptions.lik_init
DynareOptions.lik_init = 1;
Pstar = lyapunov_symm(T,R*Q*R',DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold);
end
Pinf = [];
Pinf = [];
a = zeros(mm,1);
Zflag = 0;
otherwise
@ -471,9 +471,13 @@ if analytic_derivation,
no_DLIK = 0;
full_Hess = analytic_derivation==2;
asy_Hess = analytic_derivation==-2;
outer_product_gradient = analytic_derivation==-1;
if asy_Hess,
analytic_derivation=1;
end
if outer_product_gradient,
analytic_derivation=1;
end
DLIK = [];
AHess = [];
if nargin<8 || isempty(derivatives_info)
@ -516,7 +520,7 @@ if analytic_derivation,
DT = DT(iv,iv,:);
DOm = DOm(iv,iv,:);
DYss = DYss(iv,:);
DH=zeros([size(H),length(xparam1)]);
DH=zeros([length(H),length(H),length(xparam1)]);
DQ=zeros([size(Q),length(xparam1)]);
DP=zeros([size(T),length(xparam1)]);
if full_Hess,
@ -578,7 +582,7 @@ if analytic_derivation,
end
end
if analytic_derivation==1,
analytic_deriv_info={analytic_derivation,DT,DYss,DOm,DH,DP};
analytic_deriv_info={analytic_derivation,DT,DYss,DOm,DH,DP,asy_Hess};
else
analytic_deriv_info={analytic_derivation,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P};
end
@ -614,6 +618,8 @@ if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter
if analytic_derivation,
LIK1=LIK;
LIK=LIK1{1};
lik1=lik;
lik=lik1{1};
end
if isinf(LIK)
if kalman_algo == 1
@ -645,6 +651,7 @@ if (kalman_algo==2) || (kalman_algo==4)
if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
H = diag(H);
mmm = mm;
clear tmp
if analytic_derivation,
for j=1:pp,
tmp(j,:)=DH(j,j,:);
@ -676,6 +683,8 @@ if (kalman_algo==2) || (kalman_algo==4)
if analytic_derivation,
LIK1=LIK;
LIK=LIK1{1};
lik1=lik;
lik=lik1{1};
end
if DynareOptions.lik_init==3
LIK = LIK+dLIK;
@ -690,13 +699,17 @@ if analytic_derivation
DLIK = LIK1{2};
% [DLIK] = score(T,R,Q,H,Pstar,Y,DT,DYss,DOm,DH,DP,start,Z,kalman_tol,riccati_tol);
end
if full_Hess,
if full_Hess ,
Hess = -LIK1{3};
% [Hess, DLL] = get_Hessian(T,R,Q,H,Pstar,Y,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P,start,Z,kalman_tol,riccati_tol);
% Hess0 = getHessian(Y,T,DT,D2T, R*Q*transpose(R),DOm,D2Om,Z,DYss,D2Yss);
end
if asy_Hess,
[Hess] = AHessian(T,R,Q,H,Pstar,Y,DT,DYss,DOm,DH,DP,start,Z,kalman_tol,riccati_tol);
% if ~((kalman_algo==2) || (kalman_algo==4)),
% [Hess] = AHessian(T,R,Q,H,Pstar,Y,DT,DYss,DOm,DH,DP,start,Z,kalman_tol,riccati_tol);
% else
Hess = LIK1{3};
% end
end
end
@ -725,6 +738,11 @@ if analytic_derivation
if no_DLIK==0
DLIK = DLIK - dlnprior';
end
if outer_product_gradient,
dlik = lik1{2};
dlik=[- dlnprior; dlik(start:end,:)];
Hess = dlik'*dlik;
end
else
lnprior = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4);
end

19
matlab/dynare_estimation_1.m
View File

@ -213,22 +213,29 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
else
flag = 1;
end
if ~exist('igg','var'), % by M. Ratto
hh=[];
gg=[];
igg=[];
end % by M. Ratto
if isfield(options_,'ftol')
crit = options_.ftol;
else
crit = 1.e-5;
end
if options_.analytic_derivation,
analytic_grad=1;
ana_deriv = options_.analytic_derivation;
options_.analytic_derivation = -1;
crit = 1.e-7;
flag = 0;
else
analytic_grad=0;
end
if isfield(options_,'nit')
nit = options_.nit;
else
nit=1000;
end
[xparam1,hh,gg,fval,invhess] = newrat(objective_function,xparam1,hh,gg,igg,crit,nit,flag,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
[xparam1,hh,gg,fval,invhess] = newrat(objective_function,xparam1,analytic_grad,crit,nit,flag,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
if options_.analytic_derivation,
options_.analytic_derivation = ana_deriv;
end
parameter_names = bayestopt_.name;
save([M_.fname '_mode.mat'],'xparam1','hh','gg','fval','invhess','parameter_names');
case 6

7
matlab/dynare_identification.m
View File

@ -123,10 +123,13 @@ if isempty(dataset_),
dataset_.info.varobs = options_.varobs;
dataset_.rawdata = [];
dataset_.missing.state = 0;
dataset_.missing.aindex = [];
for jdata=1:periods,
temp1{jdata}=[1:dataset_.info.nvobs]';
end
dataset_.missing.aindex = temp1;
dataset_.missing.vindex = [];
dataset_.missing.number_of_observations = [];
dataset_.missing.no_more_missing_observations = [];
dataset_.missing.no_more_missing_observations = 1;
dataset_.descriptive.mean = [];
dataset_.data = [];

23
matlab/dynare_sensitivity.m
View File

@ -79,6 +79,11 @@ if options_gsa.identification,
options_gsa = set_default_option(options_gsa,'ar',1);
options_gsa = set_default_option(options_gsa,'useautocorr',0);
options_.ar = options_gsa.ar;
if options_gsa.morris==0,
disp('The option morris = 0 is no longer supported (Type I errors)')
disp('This option is reset at morris = 2 (local identification analysis).')
options_gsa.morris=2;
end
if options_gsa.morris==2,
if isfield(options_,'options_ident'),
options_.options_ident.load_ident_files = options_gsa.load_ident_files;
@ -122,7 +127,13 @@ if options_gsa.redform,
options_gsa.ppost=0;
end
if options_gsa.morris==1 || options_gsa.morris==3,
if options_gsa.morris>2,
disp('The option morris = 3 is no longer supported')
disp('the option is reset at morris = 1 .')
options_gsa.morris=1;
end
if options_gsa.morris==1,
if ~options_gsa.identification,
options_gsa.redform=1;
end
@ -134,12 +145,12 @@ if options_gsa.morris==1 || options_gsa.morris==3,
options_gsa.load_rmse=0;
options_gsa.alpha2_stab=1;
options_gsa.ksstat=1;
if options_gsa.morris==3,
options_gsa = set_default_option(options_gsa,'Nsam',256);
OutputDirectoryName = CheckPath('gsa/identif',M_.dname);
else
% if options_gsa.morris==3,
% options_gsa = set_default_option(options_gsa,'Nsam',256);
% OutputDirectoryName = CheckPath('gsa/identif',M_.dname);
% else
OutputDirectoryName = CheckPath('gsa/screen',M_.dname);
end
% end
else
OutputDirectoryName = CheckPath('gsa',M_.dname);
end

2
matlab/getH.m
View File

@ -338,7 +338,7 @@ if kronflag==1, % kronecker products
elseif kronflag==-1, % perturbation
fun = 'thet2tau';
params0 = M_.params;
H = fjaco(fun,[sqrt(diag(M_.Sigma_e(indexo,indexo))); M_.params(indx)],indx,indexo);
H = fjaco(fun,[sqrt(diag(M_.Sigma_e(indexo,indexo))); M_.params(indx)], M_, oo_, indx, indexo);
assignin('base','M_', M_);
assignin('base','oo_', oo_);

4
matlab/getJJ.m
View File

@ -28,10 +28,10 @@ warning('off','MATLAB:divideByZero')
if kronflag == -1,
fun = 'thet2tau';
params0 = M_.params;
JJ = fjaco(fun,[sqrt(diag(M_.Sigma_e(indexo,indexo))); M_.params(indx)],indx,indexo,1,mf,nlags,useautocorr);
JJ = fjaco(fun,[sqrt(diag(M_.Sigma_e(indexo,indexo))); M_.params(indx)],M_, oo_, indx,indexo,1,mf,nlags,useautocorr);
M_.params = params0;
params0 = M_.params;
H = fjaco(fun,[sqrt(diag(M_.Sigma_e(indexo,indexo))); M_.params(indx)],indx,indexo,0,mf,nlags,useautocorr);
H = fjaco(fun,[sqrt(diag(M_.Sigma_e(indexo,indexo))); M_.params(indx)],M_, oo_, indx,indexo,0,mf,nlags,useautocorr);
M_.params = params0;
assignin('base','M_', M_);
assignin('base','oo_', oo_);

2
matlab/gsa/map_ident_.m
View File

@ -768,6 +768,8 @@ if opt_gsa.morris==1,
% eval(['print -dpdf ' OutputDirectoryName '/' fname_ '_morris_redform']);
elseif opt_gsa.morris==3,
return
np=estim_params_.np;
na=(4*np+1)*opt_gsa.Nsam;
for j=1:j0,

6
matlab/gsa/stab_map_.m
View File

@ -122,9 +122,9 @@ if fload==0,
Nsam=size(lpmat,1);
lpmat0 = lpmat(:,1:nshock);
lpmat = lpmat(:,nshock+1:end);
elseif opt_gsa.morris==3,
lpmat = prep_ide(Nsam,np,5);
Nsam=size(lpmat,1);
% elseif opt_gsa.morris==3,
% lpmat = prep_ide(Nsam,np,5);
% Nsam=size(lpmat,1);
else
if np<52 & ilptau>0,
[lpmat] = qmc_sequence(np, int64(0), 0, Nsam)';

62
matlab/gsa/thet2tau.m
View File

@ -1,62 +0,0 @@
function tau = thet2tau(params, indx, indexo, flagmoments,mf,nlags,useautocorr)
% Copyright (C) 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 <http://www.gnu.org/licenses/>.
global M_ oo_ options_
if nargin==1,
indx = [1:M_.param_nbr];
indexo = [];
end
if nargin<4,
flagmoments=0;
end
if nargin<7 | isempty(useautocorr),
useautocorr=0;
end
M_.params(indx) = params(length(indexo)+1:end);
if ~isempty(indexo)
M_.Sigma_e(indexo,indexo) = diag(params(1:length(indexo)).^2);
end
[A,B,plouf,plouf,M_,options_,oo_] = dynare_resolve(M_,options_,oo_);
if flagmoments==0,
tau = [oo_.dr.ys(oo_.dr.order_var); A(:); dyn_vech(B*M_.Sigma_e*B')];
else
GAM = lyapunov_symm(A,B*M_.Sigma_e*B',options_.qz_criterium,options_.lyapunov_complex_threshold);
k = find(abs(GAM) < 1e-12);
GAM(k) = 0;
if useautocorr,
sy = sqrt(diag(GAM));
sy = sy*sy';
sy0 = sy-diag(diag(sy))+eye(length(sy));
dum = GAM./sy0;
tau = dyn_vech(dum(mf,mf));
else
tau = dyn_vech(GAM(mf,mf));
end
for ii = 1:nlags
dum = A^(ii)*GAM;
if useautocorr,
dum = dum./sy;
end
tau = [tau;vec(dum(mf,mf))];
end
tau = [ oo_.dr.ys(oo_.dr.order_var(mf)); tau];
end

6
matlab/identification_analysis.m
View File

@ -7,7 +7,7 @@ function [ide_hess, ide_moments, ide_model, ide_lre, derivatives_info, info] = i
% o indx [array] index of estimated parameters
% o indexo [array] index of estimated shocks
% o options_ident [structure] identification options
% o data_info [structure] data info for Kalmna Filter
% o data_info [structure] data info for Kalman Filter
% o prior_exist [integer]
% =1 when prior exists and indentification is checked only for estimated params and shocks
% =0 when prior is not defined and indentification is checked for all params and shocks
@ -130,7 +130,9 @@ if info(1)==0,
options_.noprint = 1;
options_.order = 1;
options_.periods = data_info.info.ntobs+100;
options_.kalman_algo = 1;
if options_.kalman_algo > 2,
options_.kalman_algo = 1;
end
options_.analytic_derivation = -2;
info = stoch_simul(options_.varobs);
data_info.data=oo_.endo_simul(options_.varobs_id,100+1:end);

23
matlab/kalman/likelihood/computeDLIK.m
View File

@ -24,18 +24,18 @@ persistent DK DF D2K D2F
if notsteady
if Zflag
[DK,DF,DP1] = computeDKalmanZ(T,DT,DOm,P,DP,DH,Z,iF,K);
if nargout>3,
if nargout>4,
[D2K,D2F,D2P1] = computeD2KalmanZ(T,DT,D2T,D2Om,P,DP,D2P,DH,Z,iF,K,DK);
end
else
[DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,Z,iF,K);
if nargout>3,
if nargout>4,
[D2K,D2F,D2P1] = computeD2Kalman(T,DT,D2T,D2Om,P,DP,D2P,DH,Z,iF,K,DK);
end
end
else
DP1=DP;
if nargout>3,
if nargout>4,
D2P1=D2P;
end
end
@ -95,10 +95,27 @@ for ii = 1:k
end
end
end
Da(:,ii) = DT(:,:,ii)*tmp + T*dtmp(:,ii);
DLIK(ii,1) = trace( iF*DF(:,:,ii) ) + 2*Dv(:,ii)'*iF*v - v'*(iF*DF(:,:,ii)*iF)*v;
end
if nargout==4,
% Hesst(ii,jj) = getHesst_ij(v,Dv(:,ii),Dv(:,jj),0,iF,diFi,diFj,0,dFj,0);
vecDPmf = reshape(DF,[],k);
D2a = 2*Dv'*iF*Dv + (vecDPmf' * kron(iF,iF) * vecDPmf);
% for ii = 1:k
%
% diFi = -iF*DF(:,:,ii)*iF;
% for jj = 1:ii
% dFj = DF(:,:,jj);
% diFj = -iF*DF(:,:,jj)*iF;
%
% Hesst(ii,jj) = getHesst_ij(v*0,Dv(:,ii),Dv(:,jj),v*0,iF,diFi,diFj,0,-dFj,0);
% end
% end
end
% end of computeDLIK
function Hesst_ij = getHesst_ij(e,dei,dej,d2eij,iS,diSi,diSj,d2iSij,dSj,d2Sij);

52
matlab/kalman/likelihood/kalman_filter.m
View File

@ -1,4 +1,4 @@
function [LIK, likk, a, P] = kalman_filter(Y,start,last,a,P,kalman_tol,riccati_tol,presample,T,Q,R,H,Z,mm,pp,rr,Zflag,diffuse_periods,analytic_derivation,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P)
function [LIK, LIKK, a, P] = kalman_filter(Y,start,last,a,P,kalman_tol,riccati_tol,presample,T,Q,R,H,Z,mm,pp,rr,Zflag,diffuse_periods,analytic_derivation,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P)
% Computes the likelihood of a stationnary state space model.
%@info:
@ -123,14 +123,17 @@ LIK = Inf; % Default value of the log likelihood.
oldK = Inf;
notsteady = 1;
F_singular = 1;
asy_hess=0;
if analytic_derivation == 0,
DLIK=[];
Hess=[];
LIKK=[];
else
k = size(DT,3); % number of structural parameters
DLIK = zeros(k,1); % Initialization of the score.
Da = zeros(mm,k); % Derivative State vector.
dlikk = zeros(smpl,k);
if Zflag==0,
C = zeros(pp,mm);
@ -144,12 +147,17 @@ else
D2a = zeros(mm,k,k); % State vector.
d2C = zeros(pp,mm,k,k);
else
asy_hess=D2T;
Hess=[];
D2a=[];
D2T=[];
D2Yss=[];
end
if asy_hess,
Hess = zeros(k,k); % Initialization of the Hessian
end
LIK={inf,DLIK,Hess};
LIKK={likk,dlikk};
end
while notsteady && t<=last
@ -185,14 +193,15 @@ while notsteady && t<=last
if analytic_derivation==2,
[Da,DP,DLIKt,D2a,D2P, Hesst] = computeDLIK(k,tmp,Z,Zflag,v,T,K,P,iF,Da,DYss,DT,DOm,DP,DH,notsteady,D2a,D2Yss,D2T,D2Om,D2P);
else
[Da,DP,DLIKt] = computeDLIK(k,tmp,Z,Zflag,v,T,K,P,iF,Da,DYss,DT,DOm,DP,DH,notsteady);
[Da,DP,DLIKt,Hesst] = computeDLIK(k,tmp,Z,Zflag,v,T,K,P,iF,Da,DYss,DT,DOm,DP,DH,notsteady);
end
if t>presample
DLIK = DLIK + DLIKt;
if analytic_derivation==2,
if analytic_derivation==2 || asy_hess,
Hess = Hess + Hesst;
end
end
dlikk(s,:)=DLIKt;
end
a = T*tmp;
P=Ptmp;
@ -208,19 +217,31 @@ end
% Add observation's densities constants and divide by two.
likk(1:s) = .5*(likk(1:s) + pp*log(2*pi));
if analytic_derivation,
DLIK = DLIK/2;
dlikk = dlikk/2;
if analytic_derivation==2 || asy_hess,
if asy_hess==0,
Hess = Hess + tril(Hess,-1)';
end
Hess = -Hess/2;
end
end
% Call steady state Kalman filter if needed.
if t <= last
if analytic_derivation,
if analytic_derivation==2,
[tmp, likk(s+1:end)] = kalman_filter_ss(Y,t,last,a,T,K,iF,dF,Z,pp,Zflag, ...
[tmp, tmp2] = kalman_filter_ss(Y,t,last,a,T,K,iF,dF,Z,pp,Zflag, ...
analytic_derivation,Da,DT,DYss,D2a,D2T,D2Yss);
else
[tmp, likk(s+1:end)] = kalman_filter_ss(Y,t,last,a,T,K,iF,dF,Z,pp,Zflag, ...
analytic_derivation,Da,DT,DYss);
[tmp, tmp2] = kalman_filter_ss(Y,t,last,a,T,K,iF,dF,Z,pp,Zflag, ...
analytic_derivation,Da,DT,DYss,asy_hess);
end
likk(s+1:end)=tmp2{1};
dlikk(s+1:end,:)=tmp2{2};
DLIK = DLIK + tmp{2};
if analytic_derivation==2,
if analytic_derivation==2 || asy_hess,
Hess = Hess + tmp{3};
end
else
@ -229,20 +250,19 @@ if t <= last
end
% Compute minus the log-likelihood.
if presample
if presample>=diffuse_periods
likk = likk(1+(presample-diffuse_periods):end);
end
if presample>diffuse_periods,
LIK = sum(likk(1+(presample-diffuse_periods):end));
else
LIK = sum(likk);
end
LIK = sum(likk);
if analytic_derivation,
DLIK = DLIK/2;
if analytic_derivation==2,
Hess = Hess + tril(Hess,-1)';
Hess = -Hess/2;
if analytic_derivation==2 || asy_hess,
LIK={LIK, DLIK, Hess};
else
LIK={LIK, DLIK};
end
LIKK={likk, dlikk};
else
LIKK=likk;
end

28
matlab/kalman/likelihood/kalman_filter_ss.m
View File

@ -81,6 +81,7 @@ t = start; % Initialization of the time index.
likk = zeros(smpl,1); % Initialization of the vector gathering the densities.
LIK = Inf; % Default value of the log likelihood.
notsteady = 0;
asy_hess=0;
if nargin<12
analytic_derivation = 0;
end
@ -91,10 +92,16 @@ if analytic_derivation == 0,
else
k = size(DT,3); % number of structural parameters
DLIK = zeros(k,1); % Initialization of the score.
dlikk = zeros(smpl,k);
if analytic_derivation==2,
Hess = zeros(k,k); % Initialization of the Hessian
else
Hess=[];
asy_hess=D2a;
if asy_hess,
Hess = zeros(k,k); % Initialization of the Hessian
else
Hess=[];
end
end
end
@ -109,12 +116,13 @@ while t <= last
if analytic_derivation==2,
[Da,junk,DLIKt,D2a,junk2, Hesst] = computeDLIK(k,tmp,Z,Zflag,v,T,K,[],iF,Da,DYss,DT,[],[],[],notsteady,D2a,D2Yss,D2T,[],[]);
else
[Da,junk,DLIKt] = computeDLIK(k,tmp,Z,Zflag,v,T,K,[],iF,Da,DYss,DT,[],[],[],notsteady);
[Da,junk,DLIKt,Hesst] = computeDLIK(k,tmp,Z,Zflag,v,T,K,[],iF,Da,DYss,DT,[],[],[],notsteady);
end
DLIK = DLIK + DLIKt;
if analytic_derivation==2,
if analytic_derivation==2 || asy_hess,
Hess = Hess + Hesst;
end
dlikk(t-start+1,:)=DLIKt;
end
a = T*tmp;
likk(t-start+1) = transpose(v)*iF*v;
@ -129,9 +137,17 @@ likk = .5*(likk + pp*log(2*pi));
% Sum the observation's densities (minus the likelihood)
LIK = sum(likk);
if analytic_derivation==2,
LIK={LIK,DLIK,Hess};
if analytic_derivation,
dlikk = dlikk/2;
DLIK = DLIK/2;
likk = {likk, dlikk};
end
if analytic_derivation==1,
if analytic_derivation==2 || asy_hess,
if asy_hess==0,
Hess = Hess + tril(Hess,-1)';
end
Hess = -Hess/2;
LIK={LIK,DLIK,Hess};
elseif analytic_derivation==1,
LIK={LIK,DLIK};
end

21
matlab/kalman/likelihood/univariate_computeDLIK.m
View File

@ -30,7 +30,7 @@ if notsteady,
DF(j)=Z*DP(:,:,j)*Z'+DH;
DK(:,j) = (DP(:,:,j)*Z')/F-PZ*DF(j)/F^2;
end
if nargout>3
if nargout>4
D2F = zeros(k,k);
D2v = zeros(k,k);
D2K = zeros(rows(K),k,k);
@ -50,7 +50,7 @@ if notsteady,
Dv = -Da(Z,:) - DYss(Z,:);
DF = squeeze(DP(Z,Z,:))+DH';
DK = squeeze(DP(:,Z,:))/F-PZ*transpose(DF)/F^2;
if nargout>3
if nargout>4
D2v = squeeze(-D2a(Z,:,:) - D2Yss(Z,:,:));
D2F = squeeze(D2P(Z,Z,:,:));
D2K = squeeze(D2P(:,Z,:,:))/F;
@ -60,7 +60,7 @@ if notsteady,
end
end
end
if nargout>3
if nargout>4
DD2K(:,indx,:,:)=D2K;
DD2F(indx,:,:)=D2F;
end
@ -69,13 +69,13 @@ if notsteady,
else
DK = squeeze(DDK(:,indx,:));
DF = DDF(indx,:)';
if nargout>3
if nargout>4
D2K = squeeze(DD2K(:,indx,:,:));
D2F = squeeze(DD2F(indx,:,:));
end
if Zflag
Dv = -Z*Da(:,:) - Z*DYss(:,:);
if nargout>3
if nargout>4
D2v = zeros(k,k);
for j=1:k,
D2v(:,j) = -Z*D2a(:,:,j) - Z*D2Yss(:,:,j);
@ -83,7 +83,7 @@ else
end
else
Dv = -Da(Z,:) - DYss(Z,:);
if nargout>3
if nargout>4
D2v = squeeze(-D2a(Z,:,:) - D2Yss(Z,:,:));
end
end
@ -93,10 +93,15 @@ DLIK = DF/F + 2*Dv'/F*v - v^2/F^2*DF;
if nargout==6
Hesst = D2F/F-1/F^2*(DF*DF') + 2*D2v/F*v + 2*(Dv'*Dv)/F - 2*(DF*Dv)*v/F^2 ...
- v^2/F^2*D2F - 2*v/F^2*(Dv'*DF') + 2*v^2/F^3*(DF*DF');
elseif nargout==4,
D2a = 1/F^2*(DF*DF') + 2*(Dv'*Dv)/F ;
% D2a = -1/F^2*(DF*DF') + 2*(Dv'*Dv)/F + 2*v^2/F^3*(DF*DF') ...
% - 2*(DF*Dv)*v/F^2 - 2*v/F^2*(Dv'*DF');
% D2a = +2*(Dv'*Dv)/F + (DF' * DF)/F^2;
end
Da = Da + DK*v+K*Dv;
if nargout>3
if nargout>4
D2a = D2a + D2K*v;
for j=1:k,
@ -121,7 +126,7 @@ if notsteady,
DP1(:,:,j)=DP(:,:,j) - (DP(:,Z,j))*K'-PZ*DK(:,j)';
end
end
if nargout>3,
if nargout>4,
if Zflag,
for j=1:k,
D2P = D2P;

34
matlab/kalman/likelihood/univariate_kalman_filter.m
View File

@ -119,6 +119,7 @@ notsteady = 1;
oldK = Inf;
K = NaN(mm,pp);
asy_hess=0;
if analytic_derivation == 0,
DLIK=[];
@ -127,6 +128,7 @@ else
k = size(DT,3); % number of structural parameters
DLIK = zeros(k,1); % Initialization of the score.
Da = zeros(mm,k); % Derivative State vector.
dlik = zeros(smpl,k);
if Zflag==0,
C = zeros(pp,mm);
@ -140,11 +142,15 @@ else
D2a = zeros(mm,k,k); % State vector.
d2C = zeros(pp,mm,k,k);
else
asy_hess=D2T;
Hess=[];
D2a=[];
D2T=[];
D2Yss=[];
end
if asy_hess,
Hess = zeros(k,k); % Initialization of the Hessian
end
LIK={inf,DLIK,Hess};
end
@ -177,14 +183,15 @@ while notsteady && t<=last
if analytic_derivation==2,
[Da,DP,DLIKt,D2a,D2P, Hesst] = univariate_computeDLIK(k,i,z(i,:),Zflag,prediction_error,Ki,PZ,Fi,Da,DYss,DP,DH(d_index(i),:),notsteady,D2a,D2Yss,D2P);
else
[Da,DP,DLIKt] = univariate_computeDLIK(k,i,z(i,:),Zflag,prediction_error,Ki,PZ,Fi,Da,DYss,DP,DH(d_index(i),:),notsteady);
[Da,DP,DLIKt,Hesst] = univariate_computeDLIK(k,i,z(i,:),Zflag,prediction_error,Ki,PZ,Fi,Da,DYss,DP,DH(d_index(i),:),notsteady);
end
if t>presample
DLIK = DLIK + DLIKt;
if analytic_derivation==2,
if analytic_derivation==2 || asy_hess,
Hess = Hess + Hesst;
end
end
dlik(s,:)=DLIKt;
end
a = a + Ki*prediction_error;
P = P - PZ*Ki';
@ -208,19 +215,29 @@ end
% Divide by two.
lik(1:s) = .5*lik(1:s);
if analytic_derivation,
DLIK = DLIK/2;
dlik = dlik/2;
if analytic_derivation==2 || asy_hess,
% Hess = (Hess + Hess')/2;
Hess = -Hess/2;
end
end
% Call steady state univariate kalman filter if needed.
if t <= last
if analytic_derivation,
if analytic_derivation==2,
[tmp, lik(s+1:end)] = univariate_kalman_filter_ss(Y,t,last,a,P,kalman_tol,T,H,Z,pp,Zflag, ...
[tmp, tmp2] = univariate_kalman_filter_ss(Y,t,last,a,P,kalman_tol,T,H,Z,pp,Zflag, ...
analytic_derivation,Da,DT,DYss,DP,DH,D2a,D2T,D2Yss,D2P);
else
[tmp, lik(s+1:end)] = univariate_kalman_filter_ss(Y,t,last,a,P,kalman_tol,T,H,Z,pp,Zflag, ...
analytic_derivation,Da,DT,DYss,DP,DH);
[tmp, tmp2] = univariate_kalman_filter_ss(Y,t,last,a,P,kalman_tol,T,H,Z,pp,Zflag, ...
analytic_derivation,Da,DT,DYss,DP,DH,asy_hess);
end
lik(s+1:end)=tmp2{1};
dlik(s+1:end,:)=tmp2{2};
DLIK = DLIK + tmp{2};
if analytic_derivation==2,
if analytic_derivation==2 || asy_hess,
Hess = Hess + tmp{3};
end
else
@ -236,11 +253,10 @@ else
end
if analytic_derivation,
DLIK = DLIK/2;
if analytic_derivation==2,
Hess = -Hess/2;
if analytic_derivation==2 || asy_hess,
LIK={LIK, DLIK, Hess};
else
LIK={LIK, DLIK};
end
lik={lik, dlik};
end

26
matlab/kalman/likelihood/univariate_kalman_filter_ss.m
View File

@ -82,6 +82,7 @@ t = start; % Initialization of the time index.
likk = zeros(smpl,1); % Initialization of the vector gathering the densities.
LIK = Inf; % Default value of the log likelihood.
l2pi = log(2*pi);
asy_hess=0;
if nargin<12
analytic_derivation = 0;
@ -93,10 +94,16 @@ if analytic_derivation == 0,
else
k = size(DT,3); % number of structural parameters
DLIK = zeros(k,1); % Initialization of the score.
dlikk = zeros(smpl,k);
if analytic_derivation==2,
Hess = zeros(k,k); % Initialization of the Hessian
else
Hess=[];
asy_hess=D2a;
if asy_hess,
Hess = zeros(k,k); % Initialization of the Hessian
else
Hess=[];
end
end
end
@ -129,12 +136,13 @@ while t<=last
if analytic_derivation==2,
[Da,DPP,DLIKt,D2a,D2PP, Hesst] = univariate_computeDLIK(k,i,Z(i,:),Zflag,prediction_error,Ki,PPZ,Fi,Da,DYss,DPP,DH(i,:),0,D2a,D2Yss,D2PP);
else
[Da,DPP,DLIKt] = univariate_computeDLIK(k,i,Z(i,:),Zflag,prediction_error,Ki,PPZ,Fi,Da,DYss,DPP,DH(i,:),0);
[Da,DPP,DLIKt,Hesst] = univariate_computeDLIK(k,i,Z(i,:),Zflag,prediction_error,Ki,PPZ,Fi,Da,DYss,DPP,DH(i,:),0);
end
DLIK = DLIK + DLIKt;
if analytic_derivation==2,
if analytic_derivation==2 || asy_hess,
Hess = Hess + Hesst;
end
dlikk(s,:)=DLIKt;
end
end
end
@ -152,9 +160,15 @@ end
likk = .5*likk;
LIK = sum(likk);
if analytic_derivation==2,
LIK={LIK,DLIK,Hess};
if analytic_derivation,
dlikk = dlikk/2;
DLIK = DLIK/2;
likk = {likk, dlikk};
end
if analytic_derivation==1,
if analytic_derivation==2 || asy_hess,
% Hess = (Hess + Hess')/2;
Hess = -Hess/2;
LIK={LIK,DLIK,Hess};
elseif analytic_derivation==1,
LIK={LIK,DLIK};
end

4
matlab/mr_gstep.m
View File

@ -21,6 +21,10 @@ function [f0, x, ig] = mr_gstep(h1,x,func0,htol0,DynareDataset,DynareOptions,Mod
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
n=size(x,1);
if isempty(h1),
h1=DynareOptions.gradient_epsilon*ones(n,1);
end
if isempty(htol0)
htol = 1.e-6;

31
matlab/newrat.m
View File

@ -1,4 +1,4 @@
function [xparam1, hh, gg, fval, igg] = newrat(func0, x, hh, gg, igg, ftol0, nit, flagg, DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
function [xparam1, hh, gg, fval, igg] = newrat(func0, x, analytic_derivation, ftol0, nit, flagg, DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
% [xparam1, hh, gg, fval, igg] = newrat(func0, x, hh, gg, igg, ftol0, nit, flagg, varargin)
%
% Optimiser with outer product gradient and with sequences of univariate steps
@ -10,9 +10,7 @@ function [xparam1, hh, gg, fval, igg] = newrat(func0, x, hh, gg, igg, ftol0, nit
% of the log-likelihood to compute outer product gradient
%
% x = starting guess
% hh = initial Hessian [OPTIONAL]
% gg = initial gradient [OPTIONAL]
% igg = initial inverse Hessian [OPTIONAL]
% analytic_derivation = 1 if analytic derivs
% ftol0 = ending criterion for function change
% nit = maximum number of iterations
%
@ -56,14 +54,14 @@ gibbstol=length(BayesInfo.pshape)/50; %25;
% func0 = str2func([func2str(func0),'_hh']);
% func0 = func0;
fval0=feval(func0,x,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
[fval0,gg,hh]=feval(func0,x,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
fval=fval0;
% initialize mr_gstep and mr_hessian
mr_hessian(1,x,[],[],[],DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
outer_product_gradient=1;
outer_product_gradient=1;
if isempty(hh)
mr_hessian(1,x,[],[],[],DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
[dum, gg, htol0, igg, hhg, h1]=mr_hessian(0,x,func0,flagit,htol,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
if isempty(dum),
outer_product_gradient=0;
@ -85,6 +83,7 @@ else
hh0=hh;
hhg=hh;
igg=inv(hh);
h1=[];
end
H = igg;
disp(['Gradient norm ',num2str(norm(gg))])
@ -125,7 +124,11 @@ while norm(gg)>gtol && check==0 && jit<nit
if length(find(ig))<nx
ggx=ggx*0;
ggx(find(ig))=gg(find(ig));
hhx = reshape(dum,nx,nx);
if analytic_derivation,
hhx=hh;
else
hhx = reshape(dum,nx,nx);
end
iggx=eye(length(gg));
iggx(find(ig),find(ig)) = inv( hhx(find(ig),find(ig)) );
[fvala,x0,fc,retcode] = csminit1(func0,x0,fval,ggx,0,iggx,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
@ -159,6 +162,11 @@ while norm(gg)>gtol && check==0 && jit<nit
if (fval0(icount)-fval)<ftol
disp('No further improvement is possible!')
check=1;
if analytic_derivation,
[fvalx,gg,hh]=feval(func0,xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
hhg=hh;
H = inv(hh);
else
if flagit==2
hh=hh0;
elseif flagg>0
@ -173,6 +181,7 @@ while norm(gg)>gtol && check==0 && jit<nit
hh=hhg;
end
end
end
disp(['Actual dxnorm ',num2str(norm(x(:,end)-x(:,end-1)))])
disp(['FVAL ',num2str(fval)])
disp(['Improvement ',num2str(fval0(icount)-fval)])
@ -191,7 +200,7 @@ while norm(gg)>gtol && check==0 && jit<nit
disp(['Ftol ',num2str(ftol)])
disp(['Htol ',num2str(htol0)])
htol=htol_base;
if norm(x(:,icount)-xparam1)>1.e-12
if norm(x(:,icount)-xparam1)>1.e-12 && analytic_derivation==0,
try
save m1.mat x fval0 nig -append
catch
@ -225,6 +234,10 @@ while norm(gg)>gtol && check==0 && jit<nit
end
H = igg;
end
else
[fvalx,gg,hh]=feval(func0,xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
hhg=hh;
H = inv(hh);
end
disp(['Gradient norm ',num2str(norm(gg))])
ee=eig(hh);

8
matlab/thet2tau.m
View File

@ -1,4 +1,4 @@
function tau = thet2tau(params, indx, indexo, flagmoments,mf,nlags,useautocorr)
function tau = thet2tau(params, M_, oo_, indx, indexo, flagmoments,mf,nlags,useautocorr)
%
% Copyright (C) 2011 Dynare Team
@ -18,17 +18,17 @@ function tau = thet2tau(params, indx, indexo, flagmoments,mf,nlags,useautocorr)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global M_ oo_ options_
global options_
if nargin==1,
indx = [1:M_.param_nbr];
indexo = [];
end
if nargin<4,
if nargin<6,
flagmoments=0;
end
if nargin<7 || isempty(useautocorr),
if nargin<9 || isempty(useautocorr),
useautocorr=0;
end

15
tests/gsa/ls2003ide.mod
View File

@ -66,23 +66,14 @@ disp('CREATE SCREENING SAMPLE, CHECK FOR STABILITY AND PERFORM A SCREENING FOR I
disp('TYPE II ERRORS')
disp(' ')
disp('PRESS ENTER TO CONTUNUE');
pause;
pause(5);
dynare_sensitivity(identification=1, morris_nliv=6, morris_ntra=50);
disp('CREATE MC SAMPLE, CHECK FOR STABILITY AND PERFORM IDENTIFICATION ANALYSIS');
disp('TYPE I (main effects)')
disp(' ')
disp('PRESS ENTER TO CONTUNUE');
pause;
dynare_sensitivity(identification=1, morris=0);
disp('USE PREVIOUS MC SAMPLE AND PERFORM IDENTIFICATION ANALYSIS');
disp('WIth analytic derivatives')
disp(' ')
disp('PRESS ENTER TO CONTUNUE');
pause;
dynare_sensitivity(identification=1, load_stab=1, stab=0, morris=2);
pause(5);
dynare_sensitivity(identification=1, morris=2);
stoch_simul(order=1,irf=40);