Description of get

0001 function [Hess] = get_Hessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P,start,mf,kalman_tol,riccati_tol)
0002 % function [Hess] = get_Hessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P,start,mf,kalman_tol,riccati_tol)
0003 %
0004 % computes the 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 % NOTE: the derivative matrices (D2T,D2Om ...) are 4-dim. arrays with last
0011 % two dimensions equal to the number of structural parameters
0012 
0013 % Copyright (C) 2011 Dynare Team
0014 %
0015 % This file is part of Dynare.
0016 %
0017 % Dynare is free software: you can redistribute it and/or modify
0018 % it under the terms of the GNU General Public License as published by
0019 % the Free Software Foundation, either version 3 of the License, or
0020 % (at your option) any later version.
0021 %
0022 % Dynare is distributed in the hope that it will be useful,
0023 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0024 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0025 % GNU General Public License for more details.
0026 %
0027 % You should have received a copy of the GNU General Public License
0028 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0029 
0030 
0031     k = size(DT,3);                                 % number of structural parameters
0032     smpl = size(Y,2);                               % Sample size.
0033     pp   = size(Y,1);                               % Maximum number of observed variables.
0034     mm   = size(T,2);                               % Number of state variables.
0035     a    = zeros(mm,1);                             % State vector.
0036     Om   = R*Q*transpose(R);                        % Variance of R times the vector of structural innovations.
0037     t    = 0;                                       % Initialization of the time index.
0038     oldK = 0;
0039     notsteady   = 1;                                % Steady state flag.
0040     F_singular  = 1;
0041 
0042     Hess  = zeros(k,k);                             % Initialization of the Hessian
0043     Da    = zeros(mm,k);                             % State vector.
0044     Dv = zeros(length(mf),k);
0045     D2a    = zeros(mm,k,k);                             % State vector.
0046     D2v = zeros(length(mf),k,k);
0047 
0048     C = zeros(length(mf),mm);
0049     for ii=1:length(mf); C(ii,mf(ii))=1;end         % SELECTION MATRIX IN MEASUREMENT EQ. (FOR WHEN IT IS NOT CONSTANT)
0050     dC = zeros(length(mf),mm,k);
0051     d2C = zeros(length(mf),mm,k,k);
0052     
0053     s   = zeros(pp,1);                      % CONSTANT TERM IN MEASUREMENT EQ. (FOR WHEN IT IS NOT CONSTANT)
0054     ds  = zeros(pp,1,k);
0055     d2s = zeros(pp,1,k,k);
0056     
0057 %     for ii = 1:k
0058 %         DOm = DR(:,:,ii)*Q*transpose(R) + R*DQ(:,:,ii)*transpose(R) + R*Q*transpose(DR(:,:,ii));
0059 %     end
0060     
0061     while notsteady & t<smpl
0062         t  = t+1;
0063         v  = Y(:,t)-a(mf);
0064         F  = P(mf,mf) + H;
0065         if rcond(F) < kalman_tol
0066             if ~all(abs(F(:))<kalman_tol)
0067                 return
0068             else
0069                 a = T*a;
0070                 P = T*P*transpose(T)+Om;
0071             end
0072         else
0073             F_singular = 0;
0074             iF     = inv(F);
0075             K      = P(:,mf)*iF;
0076 
0077             [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K);
0078             [D2K,D2F,D2P1] = computeD2Kalman(T,DT,D2T,D2Om,P,DP,D2P,DH,mf,iF,K,DK);
0079             tmp = (a+K*v);
0080 
0081     for ii = 1:k
0082         Dv(:,ii)   = -Da(mf,ii) - DYss(mf,ii);
0083       %  dai = da(:,:,ii);
0084         dKi  = DK(:,:,ii);
0085         diFi = -iF*DF(:,:,ii)*iF;
0086         dtmpi = Da(:,ii)+dKi*v+K*Dv(:,ii);
0087             
0088         
0089         for jj = 1:ii
0090             dFj    = DF(:,:,jj);
0091             diFj   = -iF*DF(:,:,jj)*iF;
0092             dKj  = DK(:,:,jj);
0093             d2Kij  = D2K(:,:,jj,ii);
0094             d2Fij  = D2F(:,:,jj,ii);
0095             d2iFij = -diFi*dFj*iF -iF*d2Fij*iF -iF*dFj*diFi;
0096             dtmpj = Da(:,jj)+dKj*v+K*Dv(:,jj);
0097 
0098             d2vij  = -D2Yss(mf,jj,ii)  - D2a(mf,jj,ii); 
0099             d2tmpij = D2a(:,jj,ii) + d2Kij*v + dKj*Dv(:,ii) + dKi*Dv(:,jj) + K*d2vij;
0100             D2a(:,jj,ii) = D2T(:,:,jj,ii)*tmp + DT(:,:,jj)*dtmpi + DT(:,:,ii)*dtmpj + T*d2tmpij;            
0101 
0102             Hesst(ii,jj) = getHesst_ij(v,Dv(:,ii),Dv(:,jj),d2vij,iF,diFi,diFj,d2iFij,dFj,d2Fij);
0103         end
0104         Da(:,ii)   = DT(:,:,ii)*tmp + T*dtmpi;
0105     end
0106 %                     vecDPmf = reshape(DP(mf,mf,:),[],k);
0107 %                     iPmf = inv(P(mf,mf));
0108                     if t>=start
0109                         Hess = Hess + Hesst;
0110                     end
0111             a      = T*(a+K*v);                   
0112             P      = T*(P-K*P(mf,:))*transpose(T)+Om;
0113             DP     = DP1;
0114             D2P     = D2P1;
0115         end
0116         notsteady = max(max(abs(K-oldK))) > riccati_tol;
0117         oldK = K;
0118     end
0119 
0120     if F_singular
0121         error('The variance of the forecast error remains singular until the end of the sample')
0122     end
0123 
0124     
0125     if t < smpl
0126         t0 = t+1;
0127         while t < smpl
0128             t = t+1;
0129             v = Y(:,t)-a(mf);
0130             tmp = (a+K*v);
0131             for ii = 1:k,
0132                 Dv(:,ii)   = -Da(mf,ii)-DYss(mf,ii);
0133                 dKi  = DK(:,:,ii);
0134                 diFi = -iF*DF(:,:,ii)*iF;
0135                 dtmpi = Da(:,ii)+dKi*v+K*Dv(:,ii);
0136                 
0137                 for jj = 1:ii,
0138                     dFj    = DF(:,:,jj);
0139                     diFj   = -iF*DF(:,:,jj)*iF;
0140                     dKj  = DK(:,:,jj);
0141                     d2Kij  = D2K(:,:,jj,ii);
0142                     d2Fij  = D2F(:,:,jj,ii);
0143                     d2iFij = -diFi*dFj*iF -iF*d2Fij*iF -iF*dFj*diFi;
0144                     dtmpj = Da(:,jj)+dKj*v+K*Dv(:,jj);
0145                     
0146                     d2vij  = -D2Yss(mf,jj,ii)  - D2a(mf,jj,ii);
0147                     d2tmpij = D2a(:,jj,ii) + d2Kij*v + dKj*Dv(:,ii) + dKi*Dv(:,jj) + K*d2vij;
0148                     D2a(:,jj,ii) = D2T(:,:,jj,ii)*tmp + DT(:,:,jj)*dtmpi + DT(:,:,ii)*dtmpj + T*d2tmpij;            
0149                     
0150                     Hesst(ii,jj) = getHesst_ij(v,Dv(:,ii),Dv(:,jj),d2vij,iF,diFi,diFj,d2iFij,dFj,d2Fij);
0151                 end
0152                 Da(:,ii)   = DT(:,:,ii)*tmp + T*dtmpi;
0153             end
0154             if t>=start
0155                 Hess = Hess + Hesst;
0156             end
0157             a = T*(a+K*v);
0158         end
0159 %         Hess = Hess + .5*(smpl+t0-1)*(vecDPmf' * kron(iPmf,iPmf) * vecDPmf);
0160         %         for ii = 1:k;
0161         %             for jj = 1:ii
0162         %              H(ii,jj) = trace(iPmf*(.5*DP(mf,mf,ii)*iPmf*DP(mf,mf,jj) + Dv(:,ii)*Dv(:,jj)'));
0163         %             end
0164         %         end
0165     end
0166     
0167 Hess = Hess + tril(Hess,-1)';
0168 
0169 Hess = -Hess/2;  
0170 % end of main function
0171 
0172 function Hesst_ij = getHesst_ij(e,dei,dej,d2eij,iS,diSi,diSj,d2iSij,dSj,d2Sij);
0173 % computes (i,j) term in the Hessian
0174 
0175 Hesst_ij = trace(diSi*dSj + iS*d2Sij) + e'*d2iSij*e + 2*(dei'*diSj*e + dei'*iS*dej + e'*diSi*dej + e'*iS*d2eij);
0176 
0177 % end of getHesst_ij
0178 
0179 function [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K)
0180 
0181             k      = size(DT,3);
0182             tmp    = P-K*P(mf,:);
0183 
0184 for ii = 1:k
0185     DF(:,:,ii)  = DP(mf,mf,ii) + DH(:,:,ii); 
0186     DiF(:,:,ii) = -iF*DF(:,:,ii)*iF;
0187     DK(:,:,ii)  = DP(:,mf,ii)*iF + P(:,mf)*DiF(:,:,ii);
0188     Dtmp        = DP(:,:,ii) - DK(:,:,ii)*P(mf,:) - K*DP(mf,:,ii);
0189     DP1(:,:,ii) = DT(:,:,ii)*tmp*T' + T*Dtmp*T' + T*tmp*DT(:,:,ii)' + DOm(:,:,ii);
0190 end
0191 
0192 % end of computeDKalman
0193 
0194 function [d2K,d2S,d2P1] = computeD2Kalman(A,dA,d2A,d2Om,P0,dP0,d2P0,DH,mf,iF,K0,dK0);
0195 % computes the second derivatives of the Kalman matrices
0196 % note: A=T in main func.
0197         
0198             k      = size(dA,3);
0199             tmp    = P0-K0*P0(mf,:);
0200 [ns,no] = size(K0);
0201 
0202 % CPC = C*P0*C'; CPC = .5*(CPC+CPC');iF = inv(CPC);
0203 % APC = A*P0*C';
0204 % APA = A*P0*A';
0205 
0206 
0207 d2K  = zeros(ns,no,k,k);
0208 d2S  = zeros(no,no,k,k);
0209 d2P1 = zeros(ns,ns,k,k);
0210 
0211 for ii = 1:k
0212     dAi = dA(:,:,ii);
0213     dFi = dP0(mf,mf,ii);
0214     d2Omi = d2Om(:,:,ii);
0215     diFi = -iF*dFi*iF;
0216     dKi = dK0(:,:,ii);
0217     for jj = 1:k
0218         dAj = dA(:,:,jj);
0219         dFj = dP0(mf,mf,jj);
0220         d2Omj = d2Om(:,:,jj);
0221         dFj = dP0(mf,mf,jj);
0222         diFj = -iF*dFj*iF;
0223         dKj = dK0(:,:,jj);
0224 
0225         d2Aij = d2A(:,:,jj,ii);
0226         d2Pij = d2P0(:,:,jj,ii);
0227         d2Omij = d2Om(:,:,jj,ii);
0228        
0229     % second order
0230     
0231     d2Fij = d2Pij(mf,mf) ;
0232     
0233 %     d2APC = d2Aij*P0*C' + A*d2Pij*C' + A*P0*d2Cij' + dAi*dPj*C' + dAj*dPi*C' + A*dPj*dCi' + A*dPi*dCj' + dAi*P0*dCj' + dAj*P0*dCi';
0234     d2APC = d2Pij(:,mf);
0235     
0236     d2iF = -diFi*dFj*iF -iF*d2Fij*iF -iF*dFj*diFi;
0237     
0238     d2Kij= d2Pij(:,mf)*iF + P0(:,mf)*d2iF + dP0(:,mf,jj)*diFi + dP0(:,mf,ii)*diFj;
0239         
0240     d2KCP = d2Kij*P0(mf,:) + K0*d2Pij(mf,:) + dKi*dP0(mf,:,jj) + dKj*dP0(mf,:,ii) ;
0241     
0242     dtmpi        = dP0(:,:,ii) - dK0(:,:,ii)*P0(mf,:) - K0*dP0(mf,:,ii);
0243     dtmpj        = dP0(:,:,jj) - dK0(:,:,jj)*P0(mf,:) - K0*dP0(mf,:,jj);
0244     d2tmp = d2Pij - d2KCP;
0245 
0246     d2AtmpA = d2Aij*tmp*A' + A*d2tmp*A' + A*tmp*d2Aij' + dAi*dtmpj*A' + dAj*dtmpi*A' + A*dtmpj*dAi' + A*dtmpi*dAj' + dAi*tmp*dAj' + dAj*tmp*dAi';
0247 
0248     d2K(:,:,ii,jj)  = d2Kij; %#ok<NASGU>
0249     d2P1(:,:,ii,jj) = d2AtmpA  + d2Omij;  %#ok<*NASGU>
0250     d2S(:,:,ii,jj)  = d2Fij;
0251 %     d2iS(:,:,ii,jj) = d2iF;
0252     end
0253 end
0254 
0255 % end of computeD2Kalman
0256 
0257 
0258
get_Hessian

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
