Commit 2fecf994 authored by Marco Ratto's avatar Marco Ratto
Browse files

1) Extended optimizer = 5 for analytic derivatives;

2) Start adapting identification routines to allow computation of analytic asymptotic Hessian with KF routines
parent 45bc5c24
......@@ -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) );
......
......@@ -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)
......@@ -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
......@@ -676,6 +682,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 +698,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 +737,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
......
......@@ -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
......
......@@ -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 = [];
......
......@@ -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
......
......@@ -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);
......
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,6 +123,7 @@ LIK = Inf; % Default value of the log likelihood.
oldK = Inf;
notsteady = 1;
F_singular = 1;
asy_hess=0;
if analytic_derivation == 0,
DLIK=[];
......@@ -131,6 +132,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.
dlikk = zeros(smpl,k);
if Zflag==0,
C = zeros(pp,mm);
......@@ -144,12 +146,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 +192,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 +216,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 +249,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
......@@ -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
......@@ -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;
......
......@@ -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};