Description of AHessian

0001 function [AHess, DLIK, LIK] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,riccati_tol)
0002 % function [AHess, DLIK, LIK] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,riccati_tol)
0003 %
0004 % computes the asymptotic hessian matrix of the log-likelihood function of
0005 % a state space model (notation as in kalman_filter.m in DYNARE
0006 % Thanks to  Nikolai Iskrev
0007 %
0008 % NOTE: the derivative matrices (DT,DR ...) are 3-dim. arrays with last
0009 % dimension equal to the number of structural parameters
0010 
0011 % Copyright (C) 2011 Dynare Team
0012 %
0013 % This file is part of Dynare.
0014 %
0015 % Dynare is free software: you can redistribute it and/or modify
0016 % it under the terms of the GNU General Public License as published by
0017 % the Free Software Foundation, either version 3 of the License, or
0018 % (at your option) any later version.
0019 %
0020 % Dynare is distributed in the hope that it will be useful,
0021 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0022 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0023 % GNU General Public License for more details.
0024 %
0025 % You should have received a copy of the GNU General Public License
0026 % along with Dynare.  If not, see <http://www.gnu.org/licen
0027 
0028 
0029     k = size(DT,3);                                 % number of structural parameters
0030     smpl = size(Y,2);                               % Sample size.
0031     pp   = size(Y,1);                               % Maximum number of observed variables.
0032     mm   = size(T,2);                               % Number of state variables.
0033     a    = zeros(mm,1);                             % State vector.
0034     Om   = R*Q*transpose(R);                        % Variance of R times the vector of structural innovations.
0035     t    = 0;                                       % Initialization of the time index.
0036     oldK = 0;
0037     notsteady   = 1;                                % Steady state flag.
0038     F_singular  = 1;
0039 
0040 lik  = zeros(smpl,1);                           % Initialization of the vector gathering the densities.
0041 LIK  = Inf;                                     % Default value of the log likelihood.
0042 if nargout > 1,
0043     DLIK  = zeros(k,1);                             % Initialization of the score.
0044 end
0045     AHess  = zeros(k,k);                             % Initialization of the Hessian
0046     Da    = zeros(mm,k);                             % State vector.
0047     Dv = zeros(length(mf),k);
0048     
0049 %     for ii = 1:k
0050 %         DOm = DR(:,:,ii)*Q*transpose(R) + R*DQ(:,:,ii)*transpose(R) + R*Q*transpose(DR(:,:,ii));
0051 %     end
0052     
0053     while notsteady & t<smpl
0054         t  = t+1;
0055         v  = Y(:,t)-a(mf);
0056         F  = P(mf,mf) + H;
0057         if rcond(F) < kalman_tol
0058             if ~all(abs(F(:))<kalman_tol)
0059                 return
0060             else
0061                 a = T*a;
0062                 P = T*P*transpose(T)+Om;
0063             end
0064         else
0065             F_singular = 0;
0066             iF     = inv(F);
0067             K      = P(:,mf)*iF;
0068             lik(t) = log(det(F))+transpose(v)*iF*v;
0069 
0070             [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K);
0071             
0072                     for ii = 1:k
0073                         Dv(:,ii)   = -Da(mf,ii) - DYss(mf,ii);
0074                         Da(:,ii)   = DT(:,:,ii)*(a+K*v) + T*(Da(:,ii)+DK(:,:,ii)*v + K*Dv(:,ii));
0075                         if t>=start && nargout > 1
0076                             DLIK(ii,1)  = DLIK(ii,1) + trace( iF*DF(:,:,ii) ) + 2*Dv(:,ii)'*iF*v - v'*(iF*DF(:,:,ii)*iF)*v;
0077                         end
0078                     end
0079                     vecDPmf = reshape(DP(mf,mf,:),[],k);
0080 %                     iPmf = inv(P(mf,mf));
0081                     if t>=start
0082                         AHess = AHess + Dv'*iF*Dv + .5*(vecDPmf' * kron(iF,iF) * vecDPmf);
0083                     end
0084             a      = T*(a+K*v);                   
0085             P      = T*(P-K*P(mf,:))*transpose(T)+Om;
0086             DP     = DP1;
0087         end
0088         notsteady = max(max(abs(K-oldK))) > riccati_tol;
0089         oldK = K;
0090     end
0091 
0092     if F_singular
0093         error('The variance of the forecast error remains singular until the end of the sample')
0094     end
0095 
0096     
0097     if t < smpl
0098         t0 = t+1;
0099         while t < smpl
0100             t = t+1;
0101             v = Y(:,t)-a(mf);
0102                       for ii = 1:k
0103                         Dv(:,ii)   = -Da(mf,ii)-DYss(mf,ii);
0104                         Da(:,ii)   = DT(:,:,ii)*(a+K*v) + T*(Da(:,ii)+DK(:,:,ii)*v + K*Dv(:,ii));
0105                 if t>=start && nargout >1
0106                    DLIK(ii,1)  = DLIK(ii,1) + trace( iF*DF(:,:,ii) ) + 2*Dv(:,ii)'*iF*v - v'*(iF*DF(:,:,ii)*iF)*v;
0107                 end
0108                     end
0109              if t>=start
0110                 AHess = AHess + Dv'*iF*Dv; 
0111              end   
0112             a = T*(a+K*v);
0113         lik(t) = transpose(v)*iF*v;
0114         end
0115         AHess = AHess + .5*(smpl+t0-1)*(vecDPmf' * kron(iF,iF) * vecDPmf);
0116         if nargout > 1
0117         for ii = 1:k
0118 %             DLIK(ii,1)  = DLIK(ii,1) + (smpl-t0+1)*trace( iF*DF(:,:,ii) );
0119         end
0120         end
0121         lik(t0:smpl) = lik(t0:smpl) + log(det(F));
0122 %         for ii = 1:k;
0123 %             for jj = 1:ii
0124 %              H(ii,jj) = trace(iPmf*(.5*DP(mf,mf,ii)*iPmf*DP(mf,mf,jj) + Dv(:,ii)*Dv(:,jj)'));
0125 %             end
0126 %         end
0127     end    
0128     
0129 AHess = -AHess;  
0130 if nargout > 1,
0131     DLIK = DLIK/2;
0132 end
0133 % adding log-likelihhod constants
0134 lik = (lik + pp*log(2*pi))/2;
0135 
0136 LIK = sum(lik(start:end)); % Minus the log-likelihood.
0137 % end of main function
0138     
0139 function [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K)
0140 
0141             k      = size(DT,3);
0142             tmp    = P-K*P(mf,:);
0143 
0144 for ii = 1:k
0145     DF(:,:,ii)  = DP(mf,mf,ii) + DH(:,:,ii); 
0146     DiF(:,:,ii) = -iF*DF(:,:,ii)*iF;
0147     DK(:,:,ii)  = DP(:,mf,ii)*iF + P(:,mf)*DiF(:,:,ii);
0148     Dtmp        = DP(:,:,ii) - DK(:,:,ii)*P(mf,:) - K*DP(mf,:,ii);
0149     DP1(:,:,ii) = DT(:,:,ii)*tmp*T' + T*Dtmp*T' + T*tmp*DT(:,:,ii)' + DOm(:,:,ii);
0150 end
0151 
0152 % end of computeDKalman
0153 
0154 
0155
AHessian

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
