Estimation with analytic scores and hessian;

This includes re-setting the list of output arguments in objective functions
Added test function

matlab/DsgeVarLikelihood.m

@@ -1,4 +1,4 @@
-function [fval,exit_flag,info,PHI,SIGMAu,iXX,prior] = DsgeVarLikelihood(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
+function [fval,grad,hess,exit_flag,info,PHI,SIGMAu,iXX,prior] = DsgeVarLikelihood(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
 % Evaluates the posterior kernel of the bvar-dsge model.
 %
 % INPUTS
@@ -37,6 +37,9 @@ function [fval,exit_flag,info,PHI,SIGMAu,iXX,prior] = DsgeVarLikelihood(xparam1,
 % Declaration of the persistent variables.
 persistent penalty dsge_prior_weight_idx
 
+grad=[];
+hess=[];
+
 % Initialization of the penalty
 if ~nargin || isempty(penalty)
     penalty = 1e8;

matlab/csminwel1.m

@@ -63,7 +63,7 @@ snit=100;
 %   stailstr=[' P' num2str(i) stailstr];
 %end
 
-[f0,cost_flag] = feval(fcn,x0,varargin{:});
+[f0,junk1,junk2,cost_flag] = feval(fcn,x0,varargin{:});
 
 if ~cost_flag
     disp('Bad initial parameter.') 
@@ -79,8 +79,11 @@ if NumGrad
       case 5
         [g,badg] = numgrad5(fcn, f0, x0, epsilon, varargin{:});    
     end
-else
+elseif ischar(grad)
     [g,badg] = feval(grad,x0,varargin{:});
+else
+    g=junk1;
+    badg=0;
 end
 retcode3=101;
 x=x0;
@@ -126,8 +129,11 @@ while ~done
                   case 5
                     [g1,badg1] = numgrad5(fcn, f1, x1, epsilon, varargin{:});             
                 end
-            else
+            elseif ischar(grad),
                 [g1 badg1] = feval(grad,x1,varargin{:});
+            else
+                [junk1,g1,junk2, cost_flag] = feval(fcn,x1,varargin{:});
+                badg1 = ~cost_flag;
             end
             wall1=badg1;
             % g1
@@ -160,8 +166,11 @@ while ~done
                           case 5
                             [g2,badg2] = numgrad5(fcn, f2, x2, epsilon, varargin{:});                   
                         end
-                    else
+                    elseif ischar(grad),
                         [g2 badg2] = feval(grad,x2,varargin{:});
+                    else
+                        [junk1,g2,junk2, cost_flag] = feval(fcn,x1,varargin{:});
+                        badg2 = ~cost_flag;
                     end
                     wall2=badg2;
                     % g2
@@ -195,8 +204,11 @@ while ~done
                                   case 5
                                     [g3,badg3] = numgrad5(fcn, f3, x3, epsilon, varargin{:});                         
                                 end
-                            else
+                            elseif ischar(grad),
                                 [g3 badg3] = feval(grad,x3,varargin{:});
+                            else
+                                [junk1,g3,junk2, cost_flag] = feval(fcn,x1,varargin{:});
+                                badg3 = ~cost_flag;
                             end
                             wall3=badg3;
                             % g3
@@ -259,8 +271,11 @@ while ~done
                   case 5
                     [gh,badgh] = numgrad5(fcn, fh, xh, epsilon, varargin{:});
                 end
-            else
+            elseif ischar(grad),
                 [gh badgh] = feval(grad, xh,varargin{:});
+            else
+                [junk1,gh,junk2, cost_flag] = feval(fcn,x1,varargin{:});
+                badgh = ~cost_flag;
             end
         end
         badgh=1;
@@ -273,7 +288,7 @@ while ~done
     %badgh
     stuck = (abs(fh-f) < crit);
     if (~badg) && (~badgh) && (~stuck)
-        H = bfgsi1(H,gh-g,xh-x);
+        H = bfgsi(H,gh-g,xh-x);
     end
     if Verbose
         disp('----')

matlab/dsge_likelihood.m

@@ -1,4 +1,4 @@
-function [fval,exit_flag,ys,trend_coeff,info,Model,DynareOptions,BayesInfo,DynareResults,DLIK,AHess] = dsge_likelihood(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults,derivatives_info)
+function [fval,DLIK,Hess,exit_flag,ys,trend_coeff,info,Model,DynareOptions,BayesInfo,DynareResults] = dsge_likelihood(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults,derivatives_info)
 % Evaluates the posterior kernel of a dsge model.
 
 %@info:
@@ -151,7 +151,7 @@ exit_flag   = 1;
 info        = 0;
 singularity_flag = 0;
 DLIK        = [];
-AHess       = [];
+Hess       = [];
 
 if DynareOptions.estimation_dll
     [fval,exit_flag,ys,trend_coeff,info,params,H,Q] ...
@@ -167,12 +167,11 @@ if DynareOptions.estimation_dll
 end
 
 % Set flag related to analytical derivatives.
-if nargout > 9
+analytic_derivation = DynareOptions.analytic_derivation;
-    analytic_derivation=1;
+if nargout==1,
-else
     analytic_derivation=0;
 end
-
+    
 %------------------------------------------------------------------------------
 % 1. Get the structural parameters & define penalties
 %------------------------------------------------------------------------------
@@ -183,6 +182,9 @@ if ~isequal(DynareOptions.mode_compute,1) && any(xparam1<BayesInfo.lb)
     fval = penalty+sum((BayesInfo.lb(k)-xparam1(k)).^2);
     exit_flag = 0;
     info = 41;
+    if analytic_derivation,
+        DLIK=ones(length(xparam1),1);
+    end
     return
 end
 
@@ -192,6 +194,9 @@ if ~isequal(DynareOptions.mode_compute,1) && any(xparam1>BayesInfo.ub)
     fval = penalty+sum((xparam1(k)-BayesInfo.ub(k)).^2);
     exit_flag = 0;
     info = 42;
+    if analytic_derivation,
+        DLIK=ones(length(xparam1),1);
+    end
     return
 end
 
@@ -282,11 +287,17 @@ if info(1) == 1 || info(1) == 2 || info(1) == 5 || info(1) == 7 || info(1) == 22
     fval = penalty+1;
     info = info(1);
     exit_flag = 0;
+    if analytic_derivation,
+        DLIK=ones(length(xparam1),1);
+    end
     return
 elseif info(1) == 3 || info(1) == 4 || info(1)==6 ||info(1) == 19 || info(1) == 20 || info(1) == 21  || info(1) == 23
     fval = penalty+info(2);
     info = info(1);
     exit_flag = 0;
+    if analytic_derivation,
+        DLIK=ones(length(xparam1),1);
+    end
     return
 end
 
@@ -472,7 +483,11 @@ end
 
 if analytic_derivation
     no_DLIK = 0;
-    full_Hess = 0;
+    full_Hess = analytic_derivation==2;
+    asy_Hess = analytic_derivation==-2;
+    if asy_Hess,
+        analytic_derivation=1;
+    end
     DLIK = [];
     AHess = [];
     if nargin<8 || isempty(derivatives_info)
@@ -489,9 +504,9 @@ if analytic_derivation
         end
 
         if full_Hess,
-        [dum, DT, DOm, DYss, dum2, D2T, D2Om, D2Yss] = getH(A, B, Model,DynareResults,0,indparam,indexo);
+            [dum, DT, DOm, DYss, dum2, D2T, D2Om, D2Yss] = getH(A, B, Model,DynareResults,0,indparam,indexo);
         else
-        [dum, DT, DOm, DYss] = getH(A, B, Model,DynareResults,0,indparam,indexo);
+            [dum, DT, DOm, DYss] = getH(A, B, Model,DynareResults,0,indparam,indexo);
         end
     else
         DT = derivatives_info.DT;
@@ -576,6 +591,13 @@ if analytic_derivation
         end
         end
     end
+    if analytic_derivation==1,
+        analytic_deriv_info={analytic_derivation,DT,DYss,DOm,DH,DP};
+    else
+        analytic_deriv_info={analytic_derivation,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P};
+    end
+else
+    analytic_deriv_info={0};
 end
 
 %------------------------------------------------------------------------------
@@ -592,19 +614,9 @@ if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter
                                 a,Pstar, ...
                                 kalman_tol, riccati_tol, ...
                                 DynareOptions.presample, ...
-                                T,Q,R,H,Z,mm,pp,rr,Zflag,diffuse_periods);
+                                T,Q,R,H,Z,mm,pp,rr,Zflag,diffuse_periods, ...
-        end
+                                analytic_deriv_info{:}); 
-        if analytic_derivation
+                            
-            if no_DLIK==0
-                [DLIK] = score(T,R,Q,H,Pstar,Y,DT,DYss,DOm,DH,DP,start,Z,kalman_tol,riccati_tol);
-            end
-            if nargout==11
-                [AHess] = AHessian(T,R,Q,H,Pstar,Y,DT,DYss,DOm,DH,DP,start,Z,kalman_tol,riccati_tol);
-                if full_Hess,
-                    Hess = 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
-            end
         end
     else
         [LIK,lik] = missing_observations_kalman_filter(DynareDataset.missing.aindex,DynareDataset.missing.number_of_observations,DynareDataset.missing.no_more_missing_observations,Y,diffuse_periods+1,size(Y,2), ...
@@ -613,6 +625,10 @@ if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter
                                                DynareOptions.presample, ...
                                                T,Q,R,H,Z,mm,pp,rr,Zflag,diffuse_periods);
     end
+    if analytic_derivation,
+        LIK1=LIK;
+        LIK=LIK1{1};
+    end
     if isinf(LIK)
         if kalman_algo == 1
             kalman_algo = 2;
@@ -636,10 +652,19 @@ if (kalman_algo==2) || (kalman_algo==4)
         if isequal(H,0)
             H = zeros(pp,1);
             mmm = mm;
+            if analytic_derivation,
+                DH = zeros(pp,length(xparam1));
+            end
         else
             if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
                 H = diag(H);
                 mmm = mm;
+                if analytic_derivation,
+                    for j=1:pp,
+                        tmp(j,:)=DH(j,j,:);
+                    end
+                    DH=tmp;
+                end
             else
                 Z = [Z, eye(pp)];
                 T = blkdiag(T,zeros(pp));
@@ -651,6 +676,9 @@ if (kalman_algo==2) || (kalman_algo==4)
                 mmm   = mm+pp;
             end
         end
+        if analytic_derivation,
+            analytic_deriv_info{5}=DH;
+        end
     end
 
     [LIK, lik] = univariate_kalman_filter(DynareDataset.missing.aindex,DynareDataset.missing.number_of_observations,DynareDataset.missing.no_more_missing_observations,Y,diffuse_periods+1,size(Y,2), ...
@@ -658,7 +686,11 @@ if (kalman_algo==2) || (kalman_algo==4)
                                        DynareOptions.kalman_tol, ...
                                        DynareOptions.riccati_tol, ...
                                        DynareOptions.presample, ...
-                                       T,Q,R,H,Z,mmm,pp,rr,Zflag,diffuse_periods);
+                                       T,Q,R,H,Z,mmm,pp,rr,Zflag,diffuse_periods,analytic_deriv_info{:});
+    if analytic_derivation,
+        LIK1=LIK;
+        LIK=LIK1{1};
+    end
     if DynareOptions.lik_init==3
         LIK = LIK+dLIK;
         if analytic_derivation==0 && nargout==2,
@@ -667,6 +699,22 @@ if (kalman_algo==2) || (kalman_algo==4)
     end
 end
 
+if analytic_derivation
+    if no_DLIK==0
+        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,
+        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);
+    end
+end
+
+
 if isnan(LIK)
     info = 45;
     exit_flag = 0;
@@ -684,7 +732,7 @@ end
 if analytic_derivation
     if full_Hess,
         [lnprior, dlnprior, d2lnprior] = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4);
-        AHess = Hess + d2lnprior;
+        Hess = Hess - d2lnprior;
     else
         [lnprior, dlnprior] = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4);
     end
@@ -705,4 +753,4 @@ penalty = fval;
 if analytic_derivation==0 && nargout==2,
     lik=lik(start:end,:);
     DLIK=[-lnprior; lik(:)];
-end
+end

matlab/dynare_estimation.m

@@ -53,6 +53,12 @@ end
 
 M_.dname = dname;
 
+if options_.mode_compute && options_.analytic_derivation,
+    analytic_derivation0=options_.analytic_derivation;
+    options_.analytic_derivation=1;
+end
+    
+
 if nnobs > 1
     for i=1:nnobs
         options_.nobs = nobs(i);
@@ -64,6 +70,10 @@ else
     dynare_estimation_1(var_list,dname);
 end
 
+if options_.mode_compute && options_.analytic_derivation,
+    options_.analytic_derivation=analytic_derivation0;
+end
+
 if nnobs > 1 && horizon > 0
     mh_replic = options_.mh_replic;
     rawdata = read_variables(options_.datafile,options_.varobs,[],options_.xls_sheet,options_.xls_range);

matlab/dynare_estimation_1.m

@@ -167,6 +167,9 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
                                  'MaxFunEvals',100000,'TolFun',1e-8,'TolX',1e-6);
         if isfield(options_,'optim_opt')
             eval(['optim_options = optimset(optim_options,' options_.optim_opt ');']);
+        end
+        if options_.analytic_derivation,
+            optim_options = optimset(optim_options,'GradObj','on');
         end
             [xparam1,fval,exitflag,output,lamdba,grad,hessian_fmincon] = ...
                 fmincon(objective_function,xparam1,[],[],[],[],lb,ub,[],optim_options,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
@@ -177,14 +180,23 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
         if isfield(options_,'optim_opt')
             eval(['optim_options = optimset(optim_options,' options_.optim_opt ');']);
         end
+        if options_.analytic_derivation,
+            optim_options = optimset(optim_options,'GradObj','on');
+        end
         [xparam1,fval,exitflag] = fminunc(objective_function,xparam1,optim_options,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
       case 4
         H0 = 1e-4*eye(nx);
         crit = 1e-7;
         nit = 1000;
         verbose = 2;
+        if options_.analytic_derivation,
+            analytic_grad=1;
+        else
+            analytic_grad=[];
+        end
+
             [fval,xparam1,grad,hessian_csminwel,itct,fcount,retcodehat] = ...
-                csminwel1(objective_function,xparam1,H0,[],crit,nit,options_.gradient_method,options_.gradient_epsilon,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
+                csminwel1(objective_function,xparam1,H0,analytic_grad,crit,nit,options_.gradient_method,options_.gradient_epsilon,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
             disp(sprintf('Objective function at mode: %f',fval))
       case 5
         if isfield(options_,'hess')
@@ -350,8 +362,17 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
     end
     if ~isequal(options_.mode_compute,6) && ~isequal(options_.mode_compute,'prior')
         if options_.cova_compute == 1
-            hh = reshape(hessian(objective_function,xparam1, ...
+            if options_.analytic_derivation && strcmp(func2str(objective_function),'dsge_likelihood'),
-                                     options_.gstep,dataset_,options_,M_,estim_params_,bayestopt_,oo_),nx,nx);
+                ana_deriv = options_.analytic_derivation;
+                options_.analytic_derivation = 2;
+                [junk1, junk2, hh] = feval(objective_function,xparam1, ...
+                    dataset_,options_,M_,estim_params_,bayestopt_,oo_);
+                options_.analytic_derivation = ana_deriv;
+                
+            else
+                hh = reshape(hessian(objective_function,xparam1, ...
+                    options_.gstep,dataset_,options_,M_,estim_params_,bayestopt_,oo_),nx,nx);
+            end
         end
     end
     parameter_names = bayestopt_.name;
@@ -865,11 +886,14 @@ if (any(bayestopt_.pshape  >0 ) && options_.mh_replic) || ...
         if options_.load_mh_file && options_.use_mh_covariance_matrix
             invhess = compute_mh_covariance_matrix;
         end
+        ana_deriv = options_.analytic_derivation;
+        options_.analytic_derivation = 0;
         if options_.cova_compute
             feval(options_.posterior_sampling_method,objective_function,options_.proposal_distribution,xparam1,invhess,bounds,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
         else
             error('I Cannot start the MCMC because the hessian of the posterior kernel at the mode was not computed.')
         end
+        options_.analytic_derivation = ana_deriv;
     end
     if options_.mh_posterior_mode_estimation
         CutSample(M_, options_, estim_params_);

matlab/dynare_estimation_init.m

@@ -284,6 +284,13 @@ else
                                              bayestopt_.smoother_var_list);
 end;
 
+if options_.analytic_derivation,
+    if ~(exist('sylvester3','file')==2),        
+        dynareroot = strrep(which('dynare'),'dynare.m','');
+        addpath([dynareroot 'gensylv'])
+    end
+end
+
 % Test if the data file is declared.
 if isempty(options_.datafile)
     if gsa_flag

matlab/dynare_gradient.m

@@ -52,11 +52,11 @@ for i=1:m
     else
         h = H(:,i);
     end   
-    [Fh,flag] = feval(fcn, x+transpose(h), varargin{:});
+    [Fh,junk1,junk2,flag] = feval(fcn, x+transpose(h), varargin{:});
     if flag
         G(:,i) = (Fh-F)/epsilon;
     else
-        [Fh,flag] = feval(fcn, x-transpose(h), varargin{:});
+        [Fh,junk1,junk2,flag] = feval(fcn, x-transpose(h), varargin{:});
         if flag
             G(:,i) = (F-Fh)/epsilon;
         else

matlab/identification_analysis.m

@@ -131,11 +131,12 @@ if info(1)==0,
             options_.order = 1;
             options_.periods = data_info.info.ntobs+100;
             options_.kalman_algo = 1;
+            options_.analytic_derivation = -2;
             info = stoch_simul(options_.varobs);
             data_info.data=oo_.endo_simul(options_.varobs_id,100+1:end);
             %                         datax=data;
             derivatives_info.no_DLIK=1;
-            [fval,cost_flag,ys,trend_coeff,info,M_,options_,bayestopt_,oo_,DLIK,AHess] = dsge_likelihood(params',data_info,options_,M_,estim_params_,bayestopt_,oo_,derivatives_info);
+            [fval,DLIK,AHess,cost_flag,ys,trend_coeff,info,M_,options_,bayestopt_,oo_] = dsge_likelihood(params',data_info,options_,M_,estim_params_,bayestopt_,oo_,derivatives_info);
 %                 fval = DsgeLikelihood(xparam1,data_info,options_,M_,estim_params_,bayestopt_,oo_);
             AHess=-AHess;
             if min(eig(AHess))<0,

matlab/initial_estimation_checks.m

@@ -38,7 +38,7 @@ end
 [DynareResults.steady_state] = evaluate_steady_state(DynareResults.steady_state,Model,DynareOptions,DynareResults,DynareOptions.diffuse_filter==0);
 
 % Evaluate the likelihood.
-[fval,a,b,c,d] = feval(objective_function,xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
+[fval,junk1,junk2,a,b,c,d] = feval(objective_function,xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
 
 if DynareOptions.dsge_var
     info = b;

matlab/kalman/likelihood/computeDLIK.m

@@ -0,0 +1,273 @@
+function [Da,DP1,DLIK,D2a,D2P1,Hesst] = computeDLIK(k,tmp,Z,Zflag,v,T,K,P,iF,Da,DYss,DT,DOm,DP,DH,notsteady,D2a,D2Yss,D2T,D2Om,D2P)
+
+% 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/>.
+
+% AUTHOR(S) marco.ratto@jrc.ec.europa.eu
+
+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,
+        [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,
+        [D2K,D2F,D2P1] = computeD2Kalman(T,DT,D2T,D2Om,P,DP,D2P,DH,Z,iF,K,DK);
+    end
+end
+else
+    DP1=DP;
+    if nargout>3,
+        D2P1=D2P;
+    end
+end
+
+Dv=zeros(length(v),k);
+% D2v=zeros(length(v),k,k);
+for ii = 1:k
+    if Zflag
+        Dv(:,ii)   = -Z*Da(:,ii) - Z*DYss(:,ii);
+%         if nargout>4,
+%             for jj = 1:ii
+%                 D2v(:,jj,ii)  = -Z*D2Yss(:,jj,ii)  - Z*D2a(:,jj,ii);
+%                 D2v(:,ii,jj) = D2v(:,jj,ii);
+%             end
+%         end
+    else
+        Dv(:,ii)   = -Da(Z,ii) - DYss(Z,ii);
+%         if nargout>4,
+%             for jj = 1:ii
+%                 D2v(:,jj,ii)  = -D2Yss(Z,jj,ii)  - D2a(Z,jj,ii);
+%                 D2v(:,ii,jj) = D2v(:,jj,ii);
+%             end
+%         end
+    end
+end
+
+Hesst = zeros(k,k);
+DLIK=zeros(k,1);
+for ii = 1:k
+    %  dai = da(:,:,ii);
+    dKi  = DK(:,:,ii);
+    dtmp(:,ii) = Da(:,ii)+dKi*v+K*Dv(:,ii);
+    
+    if nargout>4,
+        diFi = -iF*DF(:,:,ii)*iF;
+        for jj = 1:ii
+            dFj    = DF(:,:,jj);
+            diFj   = -iF*DF(:,:,jj)*iF;
+            dKj  = DK(:,:,jj);
+            d2Kij  = D2K(:,:,jj,ii);
+            d2Fij  = D2F(:,:,jj,ii);
+            d2iFij = -diFi*dFj*iF -iF*d2Fij*iF -iF*dFj*diFi;
+%             dtmpj = Da(:,jj)+dKj*v+K*Dv(:,jj);
+            
+%             d2vij  = D2v(:,ii,jj);
+            if Zflag
+                d2vij  = -Z*D2Yss(:,jj,ii)  - Z*D2a(:,jj,ii);
+            else
+                d2vij  = -D2Yss(Z,jj,ii)  - D2a(Z,jj,ii);
+            end
+            d2tmpij = D2a(:,jj,ii) + d2Kij*v + dKj*Dv(:,ii) + dKi*Dv(:,jj) + K*d2vij;
+            D2a(:,jj,ii) = D2T(:,:,jj,ii)*tmp + DT(:,:,jj)*dtmp(:,ii) + DT(:,:,ii)*dtmp(:,jj) + T*d2tmpij;
+            D2a(:,ii,jj) = D2a(:,jj,ii);
+            
+            if nargout==6,
+                Hesst(ii,jj) = getHesst_ij(v,Dv(:,ii),Dv(:,jj),d2vij,iF,diFi,diFj,d2iFij,dFj,d2Fij);
+            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
+
+% end of computeDLIK
+
+function Hesst_ij = getHesst_ij(e,dei,dej,d2eij,iS,diSi,diSj,d2iSij,dSj,d2Sij);
+% computes (i,j) term in the Hessian
+
+Hesst_ij = trace(diSi*dSj + iS*d2Sij) + e'*d2iSij*e + 2*(dei'*diSj*e + dei'*iS*dej + e'*diSi*dej + e'*iS*d2eij);
+
+% end of getHesst_ij
+
+function [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,Z,iF,K)
+
+k      = size(DT,3);
+tmp    = P-K*P(Z,:);
+DF = zeros([size(iF),k]);
+DK = zeros([size(K),k]);
+DP1 = zeros([size(P),k]);
+
+for ii = 1:k
+    DF(:,:,ii)  = DP(Z,Z,ii) + DH(:,:,ii);
+    DiF = -iF*DF(:,:,ii)*iF;
+    DK(:,:,ii)  = DP(:,Z,ii)*iF + P(:,Z)*DiF;
+    Dtmp        = DP(:,:,ii) - DK(:,:,ii)*P(Z,:) - K*DP(Z,:,ii);
+    DP1(:,:,ii) = DT(:,:,ii)*tmp*T' + T*Dtmp*T' + T*tmp*DT(:,:,ii)' + DOm(:,:,ii);
+end
+
+% end of computeDKalman
+
+function [DK,DF,DP1] = computeDKalmanZ(T,DT,DOm,P,DP,DH,Z,iF,K)
+
+k      = size(DT,3);
+tmp    = P-K*Z*P;
+DF = zeros([size(iF),k]);
+DK = zeros([size(K),k]);
+DP1 = zeros([size(P),k]);
+
+for ii = 1:k
+    DF(:,:,ii)  = Z*DP(:,:,ii)*Z + DH(:,:,ii);
+    DiF = -iF*DF(:,:,ii)*iF;
+    DK(:,:,ii)  = DP(:,:,ii)*Z*iF + P(:,:)*Z*DiF;
+    Dtmp        = DP(:,:,ii) - DK(:,:,ii)*Z*P(:,:) - K*Z*DP(:,:,ii);
+    DP1(:,:,ii) = DT(:,:,ii)*tmp*T' + T*Dtmp*T' + T*tmp*DT(:,:,ii)' + DOm(:,:,ii);
+end
+
+% end of computeDKalmanZ
+
+function [d2K,d2S,d2P1] = computeD2Kalman(A,dA,d2A,d2Om,P0,dP0,d2P0,DH,Z,iF,K0,dK0);
+% computes the second derivatives of the Kalman matrices
+% note: A=T in main func.
+        
+            k      = size(dA,3);
+            tmp    = P0-K0*P0(Z,:);
+[ns,no] = size(K0);
+
+% CPC = C*P0*C'; CPC = .5*(CPC+CPC');iF = inv(CPC);
+% APC = A*P0*C';
+% APA = A*P0*A';
+
+
+d2K  = zeros(ns,no,k,k);
+d2S  = zeros(no,no,k,k);
+d2P1 = zeros(ns,ns,k,k);
+
+for ii = 1:k
+    dAi = dA(:,:,ii);
+    dFi = dP0(Z,Z,ii);
+    d2Omi = d2Om(:,:,ii);
+    diFi = -iF*dFi*iF;
+    dKi = dK0(:,:,ii);
+    for jj = 1:k
+        dAj = dA(:,:,jj);
+        dFj = dP0(Z,Z,jj);
+        d2Omj = d2Om(:,:,jj);
+        dFj = dP0(Z,Z,jj);
+        diFj = -iF*dFj*iF;
+        dKj = dK0(:,:,jj);
+
+        d2Aij = d2A(:,:,jj,ii);
+        d2Pij = d2P0(:,:,jj,ii);
+        d2Omij = d2Om(:,:,jj,ii);
+       
+    % second order
+    
+    d2Fij = d2Pij(Z,Z) ;
+    
+%     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';
+    d2APC = d2Pij(:,Z);
+    
+    d2iF = -diFi*dFj*iF -iF*d2Fij*iF -iF*dFj*diFi;
+    
+    d2Kij= d2Pij(:,Z)*iF + P0(:,Z)*d2iF + dP0(:,Z,jj)*diFi + dP0(:,Z,ii)*diFj;
+        
+    d2KCP = d2Kij*P0(Z,:) + K0*d2Pij(Z,:) + dKi*dP0(Z,:,jj) + dKj*dP0(Z,:,ii) ;
+    
+    dtmpi        = dP0(:,:,ii) - dK0(:,:,ii)*P0(Z,:) - K0*dP0(Z,:,ii);
+    dtmpj        = dP0(:,:,jj) - dK0(:,:,jj)*P0(Z,:) - K0*dP0(Z,:,jj);
+    d2tmp = d2Pij - d2KCP;
+
+    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';
+
+    d2K(:,:,ii,jj)  = d2Kij; %#ok<NASGU>
+    d2P1(:,:,ii,jj) = d2AtmpA  + d2Omij;  %#ok<*NASGU>
+    d2S(:,:,ii,jj)  = d2Fij;
+%     d2iS(:,:,ii,jj) = d2iF;
+    end
+end
+
+% end of computeD2Kalman
+
+function [d2K,d2S,d2P1] = computeD2KalmanZ(A,dA,d2A,d2Om,P0,dP0,d2P0,DH,Z,iF,K0,dK0);
+% computes the second derivatives of the Kalman matrices
+% note: A=T in main func.
+        
+            k      = size(dA,3);
+            tmp    = P0-K0*Z*P0(:,:);
+[ns,no] = size(K0);
+
+% CPC = C*P0*C'; CPC = .5*(CPC+CPC');iF = inv(CPC);
+% APC = A*P0*C';
+% APA = A*P0*A';
+
+
+d2K  = zeros(ns,no,k,k);
+d2S  = zeros(no,no,k,k);
+d2P1 = zeros(ns,ns,k,k);
+
+for ii = 1:k
+    dAi = dA(:,:,ii);
+    dFi = Z*dP0(:,:,ii)*Z;
+    d2Omi = d2Om(:,:,ii);
+    diFi = -iF*dFi*iF;
+    dKi = dK0(:,:,ii);
+    for jj = 1:k
+        dAj = dA(:,:,jj);
+        dFj = Z*dP0(:,:,jj)*Z;
+        d2Omj = d2Om(:,:,jj);
+        dFj = Z*dP0(:,:,jj)*Z;
+        diFj = -iF*dFj*iF;
+        dKj = dK0(:,:,jj);
+
+        d2Aij = d2A(:,:,jj,ii);
+        d2Pij = d2P0(:,:,jj,ii);
+        d2Omij = d2Om(:,:,jj,ii);
+       
+    % second order
+    
+    d2Fij = Z*d2Pij(:,:)*Z ;
+    
+%     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';
+    d2APC = d2Pij(:,:)*Z;
+    
+    d2iF = -diFi*dFj*iF -iF*d2Fij*iF -iF*dFj*diFi;
+    
+    d2Kij= d2Pij(:,:)*Z*iF + P0(:,:)*Z*d2iF + dP0(:,:,jj)*Z*diFi + dP0(:,:,ii)*Z*diFj;
+        
+    d2KCP = d2Kij*Z*P0(:,:) + K0*Z*d2Pij(:,:) + dKi*Z*dP0(:,:,jj) + dKj*Z*dP0(:,:,ii) ;
+    
+    dtmpi        = dP0(:,:,ii) - dK0(:,:,ii)*Z*P0(:,:) - K0*Z*dP0(:,:,ii);
+    dtmpj        = dP0(:,:,jj) - dK0(:,:,jj)*Z*P0(:,:) - K0*Z*dP0(:,:,jj);
+    d2tmp = d2Pij - d2KCP;
+
+    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';
+
+    d2K(:,:,ii,jj)  = d2Kij; %#ok<NASGU>
+    d2P1(:,:,ii,jj) = d2AtmpA  + d2Omij;  %#ok<*NASGU>
+    d2S(:,:,ii,jj)  = d2Fij;
+%     d2iS(:,:,ii,jj) = d2iF;
+    end
+end
+
+% end of computeD2KalmanZ
+
+

matlab/kalman/likelihood/kalman_filter.m

@@ -1,8 +1,8 @@
-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)
+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:
-%! @deftypefn {Function File} {[@var{LIK},@var{likk},@var{a},@var{P} ] =} kalman_filter (@var{Y}, @var{start}, @var{last}, @var{a}, @var{P}, @var{kalman_tol}, @var{riccati_tol},@var{presample},@var{T},@var{Q},@var{R},@var{H},@var{Z},@var{mm},@var{pp},@var{rr},@var{Zflag},@var{diffuse_periods})
+%! @deftypefn {Function File} {[@var{LIK},@var{likk},@var{a},@var{P} ] =} DsgeLikelihood (@var{Y}, @var{start}, @var{last}, @var{a}, @var{P}, @var{kalman_tol}, @var{riccati_tol},@var{presample},@var{T},@var{Q},@var{R},@var{H},@var{Z},@var{mm},@var{pp},@var{rr},@var{Zflag},@var{diffuse_periods})
 %! @anchor{kalman_filter}
 %! @sp 1
 %! Computes the likelihood of a stationary state space model, given initial condition for the states (mean and variance).
@@ -63,7 +63,7 @@ function [LIK, likk, a, P] = kalman_filter(Y,start,last,a,P,kalman_tol,riccati_t
 %! @sp 2
 %! @strong{This function is called by:}
 %! @sp 1
-%! @ref{dsge_likelihood}
+%! @ref{DsgeLikelihood}
 %! @sp 2
 %! @strong{This function calls:}
 %! @sp 1
@@ -93,13 +93,16 @@ function [LIK, likk, a, P] = kalman_filter(Y,start,last,a,P,kalman_tol,riccati_t
 % Set defaults.
 if nargin<17
     Zflag = 0;
-    diffuse_periods = 0;
 end
 
 if nargin<18
     diffuse_periods = 0;
 end
 
+if nargin<19
+    analytic_derivation = 0;
+end
+
 if isempty(Zflag)
     Zflag = 0;
 end
@@ -121,6 +124,34 @@ oldK = Inf;
 notsteady   = 1;
 F_singular  = 1;
 
+if  analytic_derivation == 0,
+    DLIK=[];
+    Hess=[];
+else
+    k = size(DT,3);                                 % number of structural parameters
+    DLIK  = zeros(k,1);                             % Initialization of the score.
+    Da    = zeros(mm,k);                            % Derivative State vector.
+    
+    if Zflag==0,
+        C = zeros(pp,mm);
+        for ii=1:pp; C(ii,Z(ii))=1;end         % SELECTION MATRIX IN MEASUREMENT EQ. (FOR WHEN IT IS NOT CONSTANT)
+    else
+        C=Z;
+    end
+    dC = zeros(pp,mm,k);   % either selection matrix or schur have zero derivatives
+    if analytic_derivation==2,
+        Hess  = zeros(k,k);                             % Initialization of the Hessian
+        D2a    = zeros(mm,k,k);                             % State vector.
+        d2C = zeros(pp,mm,k,k);
+    else
+        Hess=[];
+        D2a=[];
+        D2T=[];
+        D2Yss=[];
+    end
+    LIK={inf,DLIK,Hess};
+end
+
 while notsteady && t<=last
     s = t-start+1;
     if Zflag
@@ -144,12 +175,27 @@ while notsteady && t<=last
         likk(s) = log(dF)+transpose(v)*iF*v;
         if Zflag
             K = P*Z'*iF;
-            P = T*(P-K*Z*P)*transpose(T)+QQ;
+            Ptmp = T*(P-K*Z*P)*transpose(T)+QQ;
         else
             K = P(:,Z)*iF;
-            P = T*(P-K*P(Z,:))*transpose(T)+QQ;
+            Ptmp = T*(P-K*P(Z,:))*transpose(T)+QQ;
         end
-        a = T*(a+K*v);
+        tmp = (a+K*v);
+        if analytic_derivation,
+            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);
+            end
+            if t>presample
+                DLIK = DLIK + DLIKt;
+                if analytic_derivation==2,
+                    Hess = Hess + Hesst;
+                end
+            end
+        end
+        a = T*tmp;
+        P=Ptmp;
         notsteady = max(abs(K(:)-oldK))>riccati_tol;
         oldK = K(:);
     end
@@ -165,12 +211,38 @@ likk(1:s) = .5*(likk(1:s) + pp*log(2*pi));
 
 % Call steady state Kalman filter if needed.
 if t<last
-    [tmp, likk(s+1:end)] = kalman_filter_ss(Y,t,last,a,T,K,iF,dF,Z,pp,Zflag);
+    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, ...
+                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);
+        end
+        DLIK = DLIK + tmp{2};
+        if analytic_derivation==2,
+            Hess = Hess + tmp{3};
+        end
+    else
+        [tmp, likk(s+1:end)] = kalman_filter_ss(Y,t,last,a,T,K,iF,dF,Z,pp,Zflag);
+    end
 end
 
 % Compute minus the log-likelihood.
-if presample > diffuse_periods
+if presample
-    LIK = sum(likk(1+presample-diffuse_periods:end));
+    if presample>=diffuse_periods
-else
+        likk = likk(1+(presample-diffuse_periods):end);
-    LIK = sum(likk);
+    end
-end
+end
+LIK = sum(likk);
+
+if analytic_derivation,
+    DLIK = DLIK/2;
+    if analytic_derivation==2,
+        Hess = Hess + tril(Hess,-1)';
+        Hess = -Hess/2;
+        LIK={LIK, DLIK, Hess};
+    else
+        LIK={LIK, DLIK};
+    end
+end

matlab/kalman/likelihood/kalman_filter_ss.m

@@ -1,4 +1,4 @@
-function [LIK, likk, a] = kalman_filter_ss(Y,start,last,a,T,K,iF,dF,Z,pp,Zflag)
+function [LIK, likk, a] = kalman_filter_ss(Y,start,last,a,T,K,iF,dF,Z,pp,Zflag,analytic_derivation,Da,DT,DYss,D2a,D2T,D2Yss)
 % Computes the likelihood of a stationnary state space model (steady state kalman filter).
 
 %@info:
@@ -80,6 +80,23 @@ smpl = last-start+1;
 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;
+if nargin<12
+    analytic_derivation = 0;
+end
+
+if  analytic_derivation == 0,
+    DLIK=[];
+    Hess=[];
+else
+    k = size(DT,3);                                 % number of structural parameters
+    DLIK  = zeros(k,1);                             % Initialization of the score.
+    if analytic_derivation==2,
+        Hess  = zeros(k,k);                             % Initialization of the Hessian
+    else
+        Hess=[];
+    end
+end
 
 while t <= last
     if Zflag
@@ -87,7 +104,19 @@ while t <= last
     else
         v = Y(:,t)-a(Z);
     end
-    a = T*(a+K*v);
+    tmp = (a+K*v);
+    if analytic_derivation,
+        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);
+        end
+        DLIK = DLIK + DLIKt;
+        if analytic_derivation==2,
+            Hess = Hess + Hesst;
+        end
+    end
+    a = T*tmp;
     likk(t-start+1) = transpose(v)*iF*v;
     t = t+1;
 end
@@ -99,4 +128,10 @@ likk = likk + log(dF);
 likk = .5*(likk + pp*log(2*pi));
 
 % Sum the observation's densities (minus the likelihood)
-LIK = sum(likk);
+LIK = sum(likk);
+if analytic_derivation==2,
+    LIK={LIK,DLIK,Hess};
+end
+if analytic_derivation==1,
+    LIK={LIK,DLIK};
+end

matlab/kalman/likelihood/univariate_computeDLIK.m

@@ -0,0 +1,145 @@
+function [Da,DP1,DLIK,D2a,D2P,Hesst] = univariate_computeDLIK(k,indx,Z,Zflag,v,K,PZ,F,Da,DYss,DP,DH,notsteady,D2a,D2Yss,D2P)
+
+% 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/>.
+
+% AUTHOR(S) marco.ratto@jrc.ec.europa.eu
+
+persistent DDK DDF DD2K DD2F
+
+if notsteady,
+    if Zflag
+        Dv   = -Z*Da(:,:) - Z*DYss(:,:);
+        DF = zeros(k,1);
+        DK = zeros([rows(K),k]);
+        for j=1:k,
+            DF(j)=Z*DP(:,:,j)*Z'+DH;
+            DK(:,j) = (DP(:,:,j)*Z')/F-PZ*DF(j)/F^2;
+        end
+        if nargout>3
+            D2F = zeros(k,k);
+            D2v = zeros(k,k);
+            D2K = zeros(rows(K),k,k);
+            for j=1:k,
+                D2v(:,j)   = -Z*D2a(:,:,j) - Z*D2Yss(:,:,j);
+                for i=j:k,
+                    D2F(j,i)=Z*D2P(:,:,j,i)*Z';
+                    D2F(i,j)=D2F(j,i);
+                    D2K(:,j,i) = (D2P(:,:,j,i)*Z')/F-(DP(:,:,j)*Z')*DF(i)/F^2-(DP(:,:,i)*Z')*DF(j)/F^2- ...
+                        PZ*D2F(j,i)/F^2 + 2*PZ/F^3*DF(i)*DF(j);
+                    D2K(:,i,j) = D2K(:,j,i);
+                end
+            end
+        end
+        
+    else
+        Dv   = -Da(Z,:) - DYss(Z,:);
+        DF = squeeze(DP(Z,Z,:))+DH';
+        DK = squeeze(DP(:,Z,:))/F-PZ*transpose(DF)/F^2;
+        if nargout>3
+            D2v   = squeeze(-D2a(Z,:,:) - D2Yss(Z,:,:));
+            D2F = squeeze(D2P(Z,Z,:,:));
+            D2K = squeeze(D2P(:,Z,:,:))/F;
+            for j=1:k,
+                D2K(:,:,j) = D2K(:,:,j) -PZ*D2F(j,:)/F^2 - squeeze(DP(:,Z,:))*DF(j)/F^2 - ...
+                    squeeze(DP(:,Z,j))*DF'/F^2 + 2/F^3*PZ*DF'*DF(j);
+            end
+        end
+    end
+    if nargout>3
+        DD2K(:,indx,:,:)=D2K;
+        DD2F(indx,:,:)=D2F;
+    end
+    DDK(:,indx,:)=DK;
+    DDF(indx,:)=DF;
+else
+    DK = squeeze(DDK(:,indx,:));
+    DF = DDF(indx,:)';
+    if nargout>3
+        D2K = squeeze(DD2K(:,indx,:,:));
+        D2F = squeeze(DD2F(indx,:,:));
+    end
+    if Zflag
+        Dv   = -Z*Da(:,:) - Z*DYss(:,:);
+        if nargout>3
+            D2v = zeros(k,k);
+            for j=1:k,
+                D2v(:,j)   = -Z*D2a(:,:,j) - Z*D2Yss(:,:,j);
+            end
+        end
+    else
+        Dv   = -Da(Z,:) - DYss(Z,:);
+        if nargout>3
+            D2v   = squeeze(-D2a(Z,:,:) - D2Yss(Z,:,:));
+        end
+    end
+end
+
+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');
+end
+
+Da = Da + DK*v+K*Dv;
+if nargout>3
+    
+    D2a = D2a + D2K*v;
+    for j=1:k,
+        %         D2a(:,:,j) = D2a(:,:,j) + DK*Dv(j) + DK(:,j)*Dv + K*D2v(j,:);
+        for i=1:j,
+            D2a(:,j,i) = D2a(:,j,i) + DK(:,i)*Dv(j) + DK(:,j)*Dv(i) + K*D2v(j,i);
+            D2a(:,i,j) = D2a(:,j,i);
+        end
+    end
+    
+    
+end
+
+if notsteady,
+    DP1 = DP*0;
+    if Zflag,
+        for j=1:k,
+            DP1(:,:,j)=DP(:,:,j) - (DP(:,:,j)*Z')*K'-PZ*DK(:,j)';
+        end
+    else
+        for j=1:k,
+            DP1(:,:,j)=DP(:,:,j) - (DP(:,Z,j))*K'-PZ*DK(:,j)';
+        end
+    end
+    if nargout>3,
+        if Zflag,
+            for j=1:k,
+                D2P = D2P;
+            end
+        else
+            D2PZ = squeeze(D2P(:,Z,:,:));
+            DPZ = squeeze(DP(:,Z,:));
+            for j=1:k,
+                for i=j:k,
+                    D2P(:,:,j,i) = D2P(:,:,j,i) - D2PZ(:,j,i)*K' - DPZ(:,j)*DK(:,i)'- DPZ(:,i)*DK(:,j)' - PZ*squeeze(D2K(:,j,i))';
+                    D2P(:,:,i,j) = D2P(:,:,j,i);
+                end
+            end
+            
+        end
+    end
+else
+    DP1=DP;
+end
+
+

matlab/kalman/likelihood/univariate_computeDstate.m

@@ -0,0 +1,54 @@
+function [Da1,DP1,D2a,D2P] = univariate_computeDstate(k,a,P,T,Da,DP,DT,DOm,notsteady,D2a,D2P,D2T,D2Om)
+
+% 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/>.
+
+% AUTHOR(S) marco.ratto@jrc.ec.europa.eu
+
+
+DP1=DP*0;
+Da1=Da*0;
+for j=1:k,
+    Da1(:,j) = T*Da(:,j) + DT(:,:,j)*a;
+    if notsteady,
+        DP1(:,:,j) = T*DP(:,:,j)*T'+DT(:,:,j)*P*T'+T*P*DT(:,:,j)';
+    else
+        DP1=DP;
+    end
+end
+if notsteady,
+    DP1 = DP1 + DOm;
+end
+if nargout>2,
+    for j=1:k,
+        for i=1:j,
+            D2a(:,j,i) = DT(:,:,i)*Da(:,j) + DT(:,:,j)*Da(:,i) + T*D2a(:,j,i)+ D2T(:,:,j,i)*a;
+            D2a(:,i,j) = D2a(:,j,i);
+            if notsteady,
+                D2P(:,:,j,i) = T*D2P(:,:,j,i)*T' +DT(:,:,i)*DP(:,:,j)*T'+T*DP(:,:,j)*DT(:,:,i)' + ...
+                    DT(:,:,j)*DP(:,:,i)*T'+T*DP(:,:,i)*DT(:,:,j)' + ...
+                    DT(:,:,j)*P*DT(:,:,i)'+DT(:,:,i)*P*DT(:,:,j)'+ ...
+                    D2T(:,:,j,i)*P*T'+T*P*D2T(:,:,j,i)';
+                D2P(:,:,i,j) = D2P(:,:,j,i);
+            end
+        end
+    end
+    
+    if notsteady,
+        D2P = D2P + D2Om;
+    end
+end

matlab/kalman/likelihood/univariate_kalman_filter.m

@@ -120,6 +120,34 @@ notsteady = 1;
 oldK = Inf;
 K = NaN(mm,pp);
 
+if  analytic_derivation == 0,
+    DLIK=[];
+    Hess=[];
+else
+    k = size(DT,3);                                 % number of structural parameters
+    DLIK  = zeros(k,1);                             % Initialization of the score.
+    Da    = zeros(mm,k);                            % Derivative State vector.
+    
+    if Zflag==0,
+        C = zeros(pp,mm);
+        for ii=1:pp; C(ii,Z(ii))=1;end         % SELECTION MATRIX IN MEASUREMENT EQ. (FOR WHEN IT IS NOT CONSTANT)
+    else
+        C=Z;
+    end
+    dC = zeros(pp,mm,k);   % either selection matrix or schur have zero derivatives
+    if analytic_derivation==2,
+        Hess  = zeros(k,k);                             % Initialization of the Hessian
+        D2a    = zeros(mm,k,k);                             % State vector.
+        d2C = zeros(pp,mm,k,k);
+    else
+        Hess=[];
+        D2a=[];
+        D2T=[];
+        D2Yss=[];
+    end
+    LIK={inf,DLIK,Hess};
+end
+
 while notsteady && t<=last
     s = t-start+1;
     d_index = data_index{t};
@@ -144,9 +172,29 @@ while notsteady && t<=last
             if t>=no_more_missing_observations
                 K(:,i) = Ki;
             end
+            lik(s) = lik(s) + log(Fi) + prediction_error*prediction_error/Fi + l2pi;
+            if analytic_derivation,
+                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);
+                end
+                if t>presample
+                    DLIK = DLIK + DLIKt;
+                    if analytic_derivation==2,
+                        Hess = Hess + Hesst;
+                    end
+                end
+            end
             a = a + Ki*prediction_error;
             P = P - PZ*Ki';
-            lik(s) = lik(s) + log(Fi) + prediction_error*prediction_error/Fi + l2pi;
+        end
+    end
+    if analytic_derivation,        
+        if analytic_derivation==2,
+            [Da,DP,D2a,D2P] = univariate_computeDstate(k,a,P,T,Da,DP,DT,DOm,notsteady,D2a,D2P,D2T,D2Om);
+        else
+            [Da,DP] = univariate_computeDstate(k,a,P,T,Da,DP,DT,DOm,notsteady);
         end
     end
     a = T*a;
@@ -163,7 +211,21 @@ lik(1:s) = .5*lik(1:s);
 
 % Call steady state univariate kalman filter if needed.
 if t<last
-    [tmp, lik(s+1:end)] = univariate_kalman_filter_ss(Y,t,last,a,P,kalman_tol,T,H,Z,pp,Zflag);
+    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, ...
+                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);
+        end
+        DLIK = DLIK + tmp{2};
+        if analytic_derivation==2,
+            Hess = Hess + tmp{3};
+        end
+    else
+        [tmp, lik(s+1:end)] = univariate_kalman_filter_ss(Y,t,last,a,P,kalman_tol,T,H,Z,pp,Zflag);
+    end
 end
 
 % Compute minus the log-likelihood.
@@ -171,4 +233,14 @@ if presample > diffuse_periods
     LIK = sum(lik(1+presample-diffuse_periods:end));
 else
     LIK = sum(lik);
-end
+end
+
+if analytic_derivation,
+    DLIK = DLIK/2;
+    if analytic_derivation==2,
+        Hess = -Hess/2;
+        LIK={LIK, DLIK, Hess};
+    else
+        LIK={LIK, DLIK};
+    end
+end

matlab/kalman/likelihood/univariate_kalman_filter_ss.m

@@ -1,4 +1,4 @@
-function [LIK,likk,a] = univariate_kalman_filter_ss(Y,start,last,a,P,kalman_tol,T,H,Z,pp,Zflag)
+function [LIK,likk,a] = univariate_kalman_filter_ss(Y,start,last,a,P,kalman_tol,T,H,Z,pp,Zflag,analytic_derivation,Da,DT,DYss,DP,DH,D2a,D2T,D2Yss,D2P)
 % Computes the likelihood of a stationnary state space model (steady state univariate kalman filter).
 
 %@info:
@@ -83,10 +83,33 @@ likk = zeros(smpl,1);      % Initialization of the vector gathering the densitie
 LIK  = Inf;                % Default value of the log likelihood.
 l2pi = log(2*pi);
 
+if nargin<12
+    analytic_derivation = 0;
+end
+
+if  analytic_derivation == 0,
+    DLIK=[];
+    Hess=[];
+else
+    k = size(DT,3);                                 % number of structural parameters
+    DLIK  = zeros(k,1);                             % Initialization of the score.
+    if analytic_derivation==2,
+        Hess  = zeros(k,k);                             % Initialization of the Hessian
+    else
+        Hess=[];
+    end
+end
+
 % Steady state kalman filter.
 while t<=last
     s  = t-start+1;
     PP = P;
+    if analytic_derivation,
+        DPP = DP;
+        if analytic_derivation==2,
+            D2PP = D2P;
+        end
+    end
     for i=1:pp
         if Zflag
             prediction_error = Y(i,t) - Z(i,:)*a;
@@ -102,6 +125,24 @@ while t<=last
             a  = a + Ki*prediction_error;
             PP = PP - PPZ*Ki';
             likk(s) = likk(s) + log(Fi) + prediction_error*prediction_error/Fi + l2pi;
+            if analytic_derivation,
+                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);
+                end
+                DLIK = DLIK + DLIKt;
+                if analytic_derivation==2,
+                    Hess = Hess + Hesst;
+                end
+            end
+        end
+    end
+    if analytic_derivation,        
+        if analytic_derivation==2,
+            [Da,junk,D2a] = univariate_computeDstate(k,a,P,T,Da,DP,DT,[],0,D2a,D2P,D2T);
+        else
+            Da = univariate_computeDstate(k,a,P,T,Da,DP,DT,[],0);
         end
     end
     a = T*a;
@@ -110,4 +151,10 @@ end
 
 likk = .5*likk;
 
-LIK = sum(likk);
+LIK = sum(likk);
+if analytic_derivation==2,
+    LIK={LIK,DLIK,Hess};
+end
+if analytic_derivation==1,
+    LIK={LIK,DLIK};
+end

matlab/moment_function.m

@@ -1,4 +1,4 @@
-function [g,flag] = moment_function(xparams,sample_moments,dataset,options,parallel)
+function [g,grad,hess,flag] = moment_function(xparams,sample_moments,dataset,options,parallel)
 % Evaluates the moment function of the Simulated Moments Method (discrepancy between sample and
 % ).
 %
@@ -15,7 +15,7 @@ function [g,flag] = moment_function(xparams,sample_moments,dataset,options,paral
 % SPECIAL REQUIREMENTS
 %  The user has to provide a file where the moment conditions are defined.
 
-% Copyright (C) 2010-2012 Dynare Team
+% Copyright (C) 2010 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -37,14 +37,12 @@ persistent mainStream mainState
 persistent priorObjectiveValue
 
 flag = 1;
+grad=[];
+hess=[];
 
 if nargin<5
     if isempty(mainStream)
-        if matlab_ver_less_than('7.12')
+        mainStream = RandStream.getDefaultStream;
-            mainStream = RandStream.getDefaultStream;
-        else
-            mainStream = RandStream.getGlobalStream;
-        end
         mainState  = mainStream.State;
     else
         mainStream.State = mainState;

matlab/numgrad2.m

@@ -39,11 +39,11 @@ for i=1:n
     else
         tvecv=tvec(:,i);
     end
-    [fh,cost_flag] = feval(fcn, x+scale*transpose(tvecv), varargin{:});
+    [fh,junk1,junk2,cost_flag] = feval(fcn, x+scale*transpose(tvecv), varargin{:});
     if cost_flag
         g0 = (fh - f0) / (scale*delta);
     else
-        [fh,cost_flag] = feval(fcn, x-scale*transpose(tvecv), varargin{:});
+        [fh,junk1,junk2,cost_flag] = feval(fcn, x-scale*transpose(tvecv), varargin{:});
         if cost_flag
             g0 = (f0-fh) / (scale*delta);
         else

matlab/numgrad3.m

@@ -45,8 +45,8 @@ for i=1:n
     else
         tvecv=tvec(:,i);
     end
-    [f1,cost_flag1] = feval(fcn, x+scale*transpose(tvecv), varargin{:});
+    [f1,junk1,junk2,cost_flag1] = feval(fcn, x+scale*transpose(tvecv), varargin{:});
-    [f2,cost_flag2] = feval(fcn, x-scale*transpose(tvecv), varargin{:});
+    [f2,junk1,junk2,cost_flag2] = feval(fcn, x-scale*transpose(tvecv), varargin{:});
     if cost_flag1 && cost_flag2
         g0 = (f1 - f2) / (2*scale*delta);
     else

matlab/numgrad5.m

@@ -50,15 +50,15 @@ for i=1:n
     else
         tvecv=tvec(:,i);
     end
-    [f1,cost_flag1] = feval(fcn, x+scale*transpose(tvecv), varargin{:});
+    [f1,junk1,junk2,cost_flag1] = feval(fcn, x+scale*transpose(tvecv), varargin{:});
-    [f2,cost_flag2] = feval(fcn, x-scale*transpose(tvecv), varargin{:});
+    [f2,junk1,junk2,cost_flag2] = feval(fcn, x-scale*transpose(tvecv), varargin{:});
     if cost_flag1==0 || cost_flag2==0
         cost_flag3 = 0;
         cost_flag4 = 0;
         disp('numgrad:: I cannot use the five points formula!!')
     else
-        [f3,cost_flag3] = feval(fcn, x+2*scale*transpose(tvecv), varargin{:});
+        [f3,junk1,junk2,cost_flag3] = feval(fcn, x+2*scale*transpose(tvecv), varargin{:});
-        [f4,cost_flag4] = feval(fcn, x-2*scale*transpose(tvecv), varargin{:});
+        [f4,junk1,junk2,cost_flag4] = feval(fcn, x-2*scale*transpose(tvecv), varargin{:});
     end
     if cost_flag1 && cost_flag2 && cost_flag3 && cost_flag4% Five Points formula
         g0 = (8*(f1 - f2)+ f4-f3) / (12*scale*delta);

matlab/simplex_optimization_routine.m

@@ -151,6 +151,7 @@ else
     delta = 0.05;
 end
 DELTA = delta;
+zero_delta = delta/200;% To be used instead of delta if x(i) is zero.
 
 % Set max_no_improvements.
 if isfield(options,'max_no_improvements')
@@ -173,7 +174,7 @@ if verbose
     disp(' ')
 end
 initial_point = x;
-[initial_score,nopenalty] = feval(objective_function,x,varargin{:});
+[initial_score,junk1,junk2,nopenalty] = feval(objective_function,x,varargin{:});
 if ~nopenalty
     error('simplex_optimization_routine:: Initial condition is wrong!')
 else
@@ -509,7 +510,6 @@ function [v,fv,delta] = simplex_initialization(objective_function,point,point_sc
     v(:,1) = point;
     fv = zeros(n+1,1);
     fv(1) = point_score;
-    zero_delta = delta/200;% To be used instead of delta if x(i) is zero.
     if length(delta)==1
         delta = repmat(delta,n,1);
     end
@@ -522,7 +522,7 @@ function [v,fv,delta] = simplex_initialization(objective_function,point,point_sc
         end
         v(:,j+1) = y;
         x = y;
-        [fv(j+1),nopenalty_flag] = feval(objective_function,x,varargin{:});
+        [fv(j+1),junk1,junk2,nopenalty_flag] = feval(objective_function,x,varargin{:});
         if check_delta
             while ~nopenalty_flag
                 if y(j)~=0
@@ -538,7 +538,7 @@ function [v,fv,delta] = simplex_initialization(objective_function,point,point_sc
                 end
                 v(:,j+1) = y;
                 x = y;
-                [fv(j+1),nopenalty_flag] = feval(objective_function,x,varargin{:});
+                [fv(j+1),junk1,junk2,nopenalty_flag] = feval(objective_function,x,varargin{:});
             end
         end
     end

tests/analytic_derivatives/ls2003.mod

@@ -0,0 +1,69 @@
+var y y_s R pie dq pie_s de A y_obs pie_obs R_obs;
+varexo e_R e_q e_ys e_pies e_A;
+
+parameters psi1 psi2 psi3 rho_R tau alpha rr k rho_q rho_A rho_ys rho_pies;
+
+psi1 = 1.54;
+psi2 = 0.25;
+psi3 = 0.25;
+rho_R = 0.5;
+alpha = 0.3;
+rr = 2.51;
+k = 0.5;
+tau = 0.5;
+rho_q = 0.4;
+rho_A = 0.2;
+rho_ys = 0.9;
+rho_pies = 0.7;
+
+
+model(linear);
+y = y(+1) - (tau +alpha*(2-alpha)*(1-tau))*(R-pie(+1))-alpha*(tau +alpha*(2-alpha)*(1-tau))*dq(+1) + alpha*(2-alpha)*((1-tau)/tau)*(y_s-y_s(+1))-A(+1);
+pie = exp(-rr/400)*pie(+1)+alpha*exp(-rr/400)*dq(+1)-alpha*dq+(k/(tau+alpha*(2-alpha)*(1-tau)))*y+k*alpha*(2-alpha)*(1-tau)/(tau*(tau+alpha*(2-alpha)*(1-tau)))*y_s;
+pie = de+(1-alpha)*dq+pie_s;
+R = rho_R*R(-1)+(1-rho_R)*(psi1*pie+psi2*(y+alpha*(2-alpha)*((1-tau)/tau)*y_s)+psi3*de)+e_R;
+dq = rho_q*dq(-1)+e_q;
+y_s = rho_ys*y_s(-1)+e_ys;
+pie_s = rho_pies*pie_s(-1)+e_pies;
+A = rho_A*A(-1)+e_A;
+y_obs = y-y(-1)+A;
+pie_obs = 4*pie;
+R_obs = 4*R;
+end;
+
+shocks;
+	var e_R = 1.25^2;
+	var e_q = 2.5^2;
+	var e_A = 1.89;
+	var e_ys = 1.89;
+	var e_pies = 1.89;
+end;
+
+varobs y_obs R_obs pie_obs dq de;
+
+estimated_params;
+psi1 , gamma_pdf,1.5,0.5;
+psi2 , gamma_pdf,0.25,0.125;
+psi3 , gamma_pdf,0.25,0.125;
+rho_R ,beta_pdf,0.5,0.2;
+alpha ,beta_pdf,0.3,0.1;
+rr ,gamma_pdf,2.5,1;
+k , gamma_pdf,0.5,0.25;
+tau ,gamma_pdf,0.5,0.2;
+rho_q ,beta_pdf,0.4,0.2;
+rho_A ,beta_pdf,0.5,0.2;
+rho_ys ,beta_pdf,0.8,0.1;
+rho_pies,beta_pdf,0.7,0.15;
+stderr e_R,inv_gamma_pdf,1.2533,0.6551;
+stderr e_q,inv_gamma_pdf,2.5066,1.3103;
+stderr e_A,inv_gamma_pdf,1.2533,0.6551;
+stderr e_ys,inv_gamma_pdf,1.2533,0.6551;
+stderr e_pies,inv_gamma_pdf,1.88,0.9827;
+end;
+
+options_.analytic_derivation=1;
+estimation(datafile=data_ca1,first_obs=8,nobs=79,mh_nblocks=2,
+   prefilter=1,mh_jscale=0.5,mh_replic=0, mode_compute=1);
+
+identification;
+