Commit 9de8f751 authored by Sébastien Villemot's avatar Sébastien Villemot
Browse files

Merge branch 'master' into matio

parents 5f861cb0 db480787
......@@ -3032,10 +3032,11 @@ used).
(unless @code{nograph} is used).
@item graph_format = @var{FORMAT}
@itemx graph_format = ( @var{FORMAT}, @var{FORMAT}@dots{} )
@anchor{graph_format}
Specify the file format for graphs saved to disk. Possible values are
@code{eps} (the default), @code{pdf} and @code{fig} (the latter is not
available under Octave).
Specify the file format(s) for graphs saved to disk. Possible values are
@code{eps} (the default), @code{pdf} and @code{fig} (under Octave,
only @code{eps} is available).
@item noprint
Don't print anything. Useful for loops.
......@@ -3959,6 +3960,7 @@ Default value is @code{1}.
@xref{nodisplay}.
@item graph_format = @var{FORMAT}
@itemx graph_format = ( @var{FORMAT}, @var{FORMAT}@dots{} )
@xref{graph_format}.
@item lik_init = @var{INTEGER}
......@@ -4775,6 +4777,7 @@ interval. Default: @code{0.90}
@xref{nodisplay}.
@item graph_format = @var{FORMAT}
@itemx graph_format = ( @var{FORMAT}, @var{FORMAT}@dots{} )
@xref{graph_format}.
@end table
......@@ -5711,6 +5714,7 @@ Critical value for correlation @math{\rho}: plot couples of parmaters with
@xref{nodisplay}.
@item graph_format = @var{FORMAT}
@itemx graph_format = ( @var{FORMAT}, @var{FORMAT}@dots{} )
@xref{graph_format}.
@item conf_sig = @var{DOUBLE}
......@@ -5832,6 +5836,7 @@ Specify the parameter set to use. Default: @code{prior_mean}
@xref{nodisplay}.
@item graph_format = @var{FORMAT}
@itemx graph_format = ( @var{FORMAT}, @var{FORMAT}@dots{} )
@xref{graph_format}.
@end table
......
......@@ -22,7 +22,7 @@ AC_REQUIRE([AX_MATLAB_ARCH])
AC_REQUIRE([AX_MATLAB_VERSION])
AC_REQUIRE([AC_PROG_SED])
AX_COMPARE_VERSION([$MATLAB_VERSION], [lt], [6.5], [AC_MSG_ERROR([Your MATLAB is too old, please upgrade to 6.5 (R13) at least.])])
AX_COMPARE_VERSION([$MATLAB_VERSION], [lt], [7.0], [AC_MSG_ERROR([Your MATLAB is too old, please upgrade to 7.0 (R14) at least.])])
AC_MSG_CHECKING([for options to compile MEX for MATLAB])
......
......@@ -32,12 +32,16 @@ DirectoryName = [ dname '/' type ];
if ~isdir(dname)
% Make sure there isn't a file with the same name, see trac ticket #47
delete(dname)
if isfile(dname)
delete(dname)
end
mkdir('.', dname);
end
if ~isdir(DirectoryName)
% Make sure there isn't a file with the same name, see trac ticket #47
delete(DirectoryName)
if isfile(DirectoryName)
delete(DirectoryName)
end
mkdir('.',DirectoryName);
end
......@@ -34,11 +34,10 @@ function [fval,grad,hess,exit_flag,info,PHI,SIGMAu,iXX,prior] = DsgeVarLikelihoo
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global objective_function_penalty_base
% Declaration of the persistent variables.
persistent dsge_prior_weight_idx
penalty = BayesInfo.penalty;
grad=[];
hess=[];
exit_flag = [];
......@@ -85,7 +84,7 @@ exit_flag = 1;
% Return, with endogenous penalty, if some dsge-parameters are smaller than the lower bound of the prior domain.
if DynareOptions.mode_compute ~= 1 && any(xparam1 < BayesInfo.lb)
k = find(xparam1 < BayesInfo.lb);
fval = penalty+sum((BayesInfo.lb(k)-xparam1(k)).^2);
fval = objective_function_penalty_base+sum((BayesInfo.lb(k)-xparam1(k)).^2);
exit_flag = 0;
info = 41;
return;
......@@ -94,7 +93,7 @@ end
% Return, with endogenous penalty, if some dsge-parameters are greater than the upper bound of the prior domain.
if DynareOptions.mode_compute ~= 1 && any(xparam1 > BayesInfo.ub)
k = find(xparam1 > BayesInfo.ub);
fval = penalty+sum((xparam1(k)-BayesInfo.ub(k)).^2);
fval = objective_function_penalty_base+sum((xparam1(k)-BayesInfo.ub(k)).^2);
exit_flag = 0;
info = 42;
return;
......@@ -131,7 +130,7 @@ dsge_prior_weight = Model.params(dsge_prior_weight_idx);
% Is the dsge prior proper?
if dsge_prior_weight<(NumberOfParameters+NumberOfObservedVariables)/DynareDataset.info.ntobs;
fval = penalty+abs(DynareDataset.info.ntobs*dsge_prior_weight-(NumberOfParameters+NumberOfObservedVariables));
fval = objective_function_penalty_base+abs(DynareDataset.info.ntobs*dsge_prior_weight-(NumberOfParameters+NumberOfObservedVariables));
exit_flag = 0;
info = 51;
return
......@@ -148,12 +147,12 @@ end
% Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol).
if info(1) == 1 || info(1) == 2 || info(1) == 5 || info(1) == 7 || info(1) ...
== 8 || info(1) == 22 || info(1) == 24
fval = penalty+1;
fval = objective_function_penalty_base+1;
info = info(1);
exit_flag = 0;
return
elseif info(1) == 3 || info(1) == 4 || info(1) == 19 || info(1) == 20 || info(1) == 21
fval = penalty+info(2);
fval = objective_function_penalty_base+info(2);
info = info(1);
exit_flag = 0;
return
......@@ -228,7 +227,7 @@ if ~isinf(dsge_prior_weight)% Evaluation of the likelihood of the dsge-var model
if ~ispd(SIGMAu)
v = diag(SIGMAu);
k = find(v<0);
fval = penalty + sum(v(k).^2);
fval = objective_function_penalty_base + sum(v(k).^2);
info = 52;
exit_flag = 0;
return;
......
......@@ -86,7 +86,7 @@ end;
n_explod = nnz(abs(eigenvalues_) > options.qz_criterium);
result = 0;
if (nyf== n_explod) && (dr.rank == nyf)
if (nyf== n_explod) && (dr.full_rank)
result = 1;
end
......@@ -99,7 +99,7 @@ if options.noprint == 0
disp(sprintf('\nThere are %d eigenvalue(s) larger than 1 in modulus ', n_explod));
disp(sprintf('for %d forward-looking variable(s)',nyf));
disp(' ')
if dr.rank == nyf && nyf == n_explod
if result
disp('The rank condition is verified.')
else
disp('The rank conditions ISN''T verified!')
......
......@@ -653,7 +653,7 @@ else % flgresume
% load(opts.SaveFilename, 'startseed');
% randn('state', startseed);
% disp(['SEED RELOADED FROM ' opts.SaveFilename]);
startseed = randn('state'); % for retrieving in saved variables
% startseed = randn('state'); % for retrieving in saved variables
% Initialize further constants
chiN=N^0.5*(1-1/(4*N)+1/(21*N^2)); % expectation of
......
......@@ -40,7 +40,7 @@ function [fh,xh,gh,H,itct,fcount,retcodeh] = csminwel1(fcn,x0,H0,grad,crit,nit,m
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global bayestopt_
global objective_function_penalty_base
fh = [];
xh = [];
......@@ -91,7 +91,9 @@ f=f0;
H=H0;
cliff=0;
while ~done
bayestopt_.penalty = f;
% penalty for dsge_likelihood and DsgeVarLikelihood
objective_function_penalty_base = f;
g1=[]; g2=[]; g3=[];
%addition fj. 7/6/94 for control
disp('-----------------')
......
......@@ -72,7 +72,7 @@ else
chck = 0;
end;
mexErrCheck('bytecode', chck);
dr.rank = 0;
dr.full_rank = 1;
dr.eigval = [];
dr.nstatic = 0;
dr.nfwrd = 0;
......@@ -151,7 +151,6 @@ for i = 1:Size;
data(i).rank = 0;
end
dr.eigval = [dr.eigval ; data(i).eigval];
dr.rank = dr.rank + data(i).rank;
%First order approximation
if task ~= 1
[tmp1, tmp2, indx_c] = find(M_.block_structure.block(i).lead_lag_incidence(2,:));
......@@ -221,12 +220,14 @@ for i = 1:Size;
indx_c = M_.block_structure.block(i).lead_lag_incidence(3,indx_r);
data(i).eigval = 1 ./ diag(jacob(indx_r, indx_c));
data(i).rank = sum(abs(data(i).eigval) > 0);
full_rank = (rcond(jacob(indx_r, indx_c)) > 1e-9);
else
data(i).eigval = [];
data(i).rank = 0;
full_rank = 1;
end
dr.eigval = [dr.eigval ; data(i).eigval];
dr.rank = dr.rank + data(i).rank;
dr.full_rank = dr.full_rank && full_rank;
%First order approximation
if task ~= 1
if (maximum_lag > 0)
......@@ -310,11 +311,13 @@ for i = 1:Size;
if maximum_lead > 0 && n_fwrd > 0
data(i).eigval = - jacob(1 , n_pred + n - n_fwrd + 1 : n_pred + n) / jacob(1 , n_pred + n + 1 : n_pred + n + n_fwrd) ;
data(i).rank = sum(abs(data(i).eigval) > 0);
full_rank = (abs(jacob(1,n_pred+n+1: n_pred_n+n_fwrd)) > 1e-9);
else
data(i).eigval = [];
data(i).rank = 0;
full_rank = 1;
end;
dr.rank = dr.rank + data(i).rank;
dr.full_rank = dr.full_rank && full_rank;
dr.eigval = [dr.eigval ; data(i).eigval];
case 6
%% ------------------------------------------------------------------
......@@ -323,11 +326,15 @@ for i = 1:Size;
data(i).eigval = eig(- jacob(: , 1 : n_pred) / ...
jacob(: , (n_pred + n_static + 1 : n_pred + n_static + n_pred )));
data(i).rank = 0;
full_rank = (rcond(jacob(: , (n_pred + n_static + 1 : n_pred ...
+ n_static + n_pred ))) > 1e-9);
else
data(i).eigval = [];
data(i).rank = 0;
full_rank = 1;
end;
dr.eigval = [dr.eigval ; data(i).eigval];
dr.full_rank = dr.full_rank && full_rank;
if task ~= 1
if (maximum_lag > 0)
ghx = - jacob(: , 1 : n_pred) / jacob(: , n_pred + n_static + 1 : n_pred + n_static + n_pred + n_both);
......@@ -390,11 +397,14 @@ for i = 1:Size;
data(i).eigval = eig(- jacob(: , n_pred + n - n_fwrd + 1: n_pred + n))/ ...
jacob(: , n_pred + n + 1 : n_pred + n + n_fwrd);
data(i).rank = sum(abs(data(i).eigval) > 0);
full_rank = (rcond(jacob(: , n_pred + n + 1 : n_pred + n + ...
n_fwrd)) > 1e-9);
else
data(i).eigval = [];
data(i).rank = 0;
full_rank = 1;
end;
dr.rank = dr.rank + data(i).rank;
dr.full_rank = dr.full_rank && full_rank;
dr.eigval = [dr.eigval ; data(i).eigval];
case {5,8}
%% ------------------------------------------------------------------
......@@ -450,10 +460,8 @@ for i = 1:Size;
nba = nd-sdim;
if task == 1
data(i).rank = rank(w(nd-nyf+1:end,nd-nyf+1:end));
dr.rank = dr.rank + data(i).rank;
if ~exist('OCTAVE_VERSION','builtin')
data(i).eigval = eig(E,D);
end
dr.full_rank = dr.full_rank && (rcond(w(nd-nyf+1:end,nd- ...
nyf+1:end)) > 1e-9);
dr.eigval = [dr.eigval ; data(i).eigval];
end
if (verbose)
......
......@@ -130,7 +130,7 @@ function [fval,DLIK,Hess,exit_flag,ys,trend_coeff,info,Model,DynareOptions,Bayes
% AUTHOR(S) stephane DOT adjemian AT univ DASH lemans DOT FR
penalty = BayesInfo.penalty;
global objective_function_penalty_base
% Initialization of the returned variables and others...
fval = [];
......@@ -177,7 +177,7 @@ end
% Return, with endogenous penalty, if some parameters are smaller than the lower bound of the prior domain.
if ~isequal(DynareOptions.mode_compute,1) && any(xparam1<BayesInfo.lb)
k = find(xparam1<BayesInfo.lb);
fval = penalty+sum((BayesInfo.lb(k)-xparam1(k)).^2);
fval = objective_function_penalty_base+sum((BayesInfo.lb(k)-xparam1(k)).^2);
exit_flag = 0;
info = 41;
if analytic_derivation,
......@@ -189,7 +189,7 @@ end
% Return, with endogenous penalty, if some parameters are greater than the upper bound of the prior domain.
if ~isequal(DynareOptions.mode_compute,1) && any(xparam1>BayesInfo.ub)
k = find(xparam1>BayesInfo.ub);
fval = penalty+sum((xparam1(k)-BayesInfo.ub(k)).^2);
fval = objective_function_penalty_base+sum((xparam1(k)-BayesInfo.ub(k)).^2);
exit_flag = 0;
info = 42;
if analytic_derivation,
......@@ -213,7 +213,7 @@ if EstimatedParameters.ncx
a = diag(eig(Q));
k = find(a < 0);
if k > 0
fval = penalty+sum(-a(k));
fval = objective_function_penalty_base+sum(-a(k));
exit_flag = 0;
info = 43;
return
......@@ -230,7 +230,7 @@ if EstimatedParameters.ncn
a = diag(eig(H));
k = find(a < 0);
if k > 0
fval = penalty+sum(-a(k));
fval = objective_function_penalty_base+sum(-a(k));
exit_flag = 0;
info = 44;
return
......@@ -249,7 +249,7 @@ end
% Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol).
if info(1) == 1 || info(1) == 2 || info(1) == 5 || info(1) == 7 || info(1) ...
== 8 || info(1) == 22 || info(1) == 24
fval = penalty+1;
fval = objective_function_penalty_base+1;
info = info(1);
exit_flag = 0;
if analytic_derivation,
......@@ -257,7 +257,7 @@ if info(1) == 1 || info(1) == 2 || info(1) == 5 || info(1) == 7 || info(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);
fval = objective_function_penalty_base+info(2);
info = info(1);
exit_flag = 0;
if analytic_derivation,
......
function [fval,cost_flag,ys,trend_coeff,info] = dsge_posterior_kernel(xparam1,gend,data,data_index,number_of_observations,no_more_missing_observations)
% function [fval,cost_flag,ys,trend_coeff,info] = dsge_posterior_kernel(xparam1,gend,data,data_index,number_of_observations,no_more_missing_observations)
% Evaluates the posterior kernel of a dsge model.
%
% INPUTS
% xparam1 [double] vector of model parameters.
% gend [integer] scalar specifying the number of observations.
% data [double] matrix of data
% data_index [cell] cell of column vectors
% number_of_observations [integer]
% no_more_missing_observations [integer]
% OUTPUTS
% fval : value of the posterior kernel at xparam1.
% cost_flag : zero if the function returns a penalty, one otherwise.
% ys : steady state of original endogenous variables
% trend_coeff :
% info : vector of informations about the penalty:
% 41: one (many) parameter(s) do(es) not satisfied the lower bound
% 42: one (many) parameter(s) do(es) not satisfied the upper bound
%
% SPECIAL REQUIREMENTS
%
% Copyright (C) 2004-2011 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 bayestopt_ estim_params_ options_ trend_coeff_ M_ oo_
fval = [];
ys = [];
trend_coeff = [];
cost_flag = 1;
nobs = size(options_.varobs,1);
%------------------------------------------------------------------------------
% 1. Get the structural parameters & define penalties
%------------------------------------------------------------------------------
if ~isequal(options_.mode_compute,1) && any(xparam1 < bayestopt_.lb)
k = find(xparam1 < bayestopt_.lb);
fval = bayestopt_.penalty+sum((bayestopt_.lb(k)-xparam1(k)).^2);
cost_flag = 0;
info = 41;
return;
end
if ~isequal(options_.mode_compute,1) && any(xparam1 > bayestopt_.ub)
k = find(xparam1 > bayestopt_.ub);
fval = bayestopt_.penalty+sum((xparam1(k)-bayestopt_.ub(k)).^2);
cost_flag = 0;
info = 42;
return;
end
Q = M_.Sigma_e;
H = M_.H;
for i=1:estim_params_.nvx
k =estim_params_.var_exo(i,1);
Q(k,k) = xparam1(i)*xparam1(i);
end
offset = estim_params_.nvx;
if estim_params_.nvn
for i=1:estim_params_.nvn
k = estim_params_.var_endo(i,1);
H(k,k) = xparam1(i+offset)*xparam1(i+offset);
end
offset = offset+estim_params_.nvn;
end
if estim_params_.ncx
for i=1:estim_params_.ncx
k1 =estim_params_.corrx(i,1);
k2 =estim_params_.corrx(i,2);
Q(k1,k2) = xparam1(i+offset)*sqrt(Q(k1,k1)*Q(k2,k2));
Q(k2,k1) = Q(k1,k2);
end
[CholQ,testQ] = chol(Q);
if testQ %% The variance-covariance matrix of the structural innovations is not definite positive.
%% We have to compute the eigenvalues of this matrix in order to build the penalty.
a = diag(eig(Q));
k = find(a < 0);
if k > 0
fval = bayestopt_.penalty+sum(-a(k));
cost_flag = 0;
info = 43;
return
end
end
offset = offset+estim_params_.ncx;
end
if estim_params_.ncn
for i=1:estim_params_.ncn
k1 = options_.lgyidx2varobs(estim_params_.corrn(i,1));
k2 = options_.lgyidx2varobs(estim_params_.corrn(i,2));
H(k1,k2) = xparam1(i+offset)*sqrt(H(k1,k1)*H(k2,k2));
H(k2,k1) = H(k1,k2);
end
[CholH,testH] = chol(H);
if testH
a = diag(eig(H));
k = find(a < 0);
if k > 0
fval = bayestopt_.penalty+sum(-a(k));
cost_flag = 0;
info = 44;
return
end
end
offset = offset+estim_params_.ncn;
end
if estim_params_.np > 0
M_.params(estim_params_.param_vals(:,1)) = xparam1(offset+1:end);
end
M_.Sigma_e = Q;
M_.H = H;
%------------------------------------------------------------------------------
% 2. call model setup & reduction program
%------------------------------------------------------------------------------
[T,R,SteadyState,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_);
if info(1) == 1 || info(1) == 2 || info(1) == 5
fval = bayestopt_.penalty+1;
cost_flag = 0;
return
elseif info(1) == 3 || info(1) == 4 || info(1) == 20
fval = bayestopt_.penalty+info(2);%^2; % penalty power raised in DR1.m and resol already. GP July'08
cost_flag = 0;
return
end
bayestopt_.mf = bayestopt_.mf1;
if ~options_.noconstant
if options_.loglinear == 1
constant = log(SteadyState(bayestopt_.mfys));
else
constant = SteadyState(bayestopt_.mfys);
end
else
constant = zeros(nobs,1);
end
if bayestopt_.with_trend == 1
trend_coeff = zeros(nobs,1);
t = options_.trend_coeffs;
for i=1:length(t)
if ~isempty(t{i})
trend_coeff(i) = evalin('base',t{i});
end
end
trend = repmat(constant,1,gend)+trend_coeff*[1:gend];
else
trend = repmat(constant,1,gend);
end
start = options_.presample+1;
np = size(T,1);
mf = bayestopt_.mf;
no_missing_data_flag = (number_of_observations==gend*nobs);
%------------------------------------------------------------------------------
% 3. Initial condition of the Kalman filter
%------------------------------------------------------------------------------
kalman_algo = options_.kalman_algo;
if options_.lik_init == 1 % Kalman filter
if kalman_algo ~= 2
kalman_algo = 1;
end
Pstar = lyapunov_symm(T,R*Q*R',options_.qz_criterium,options_.lyapunov_complex_threshold);
Pinf = [];
elseif options_.lik_init == 2 % Old Diffuse Kalman filter
if kalman_algo ~= 2
kalman_algo = 1;
end
Pstar = 10*eye(np);
Pinf = [];
elseif options_.lik_init == 3 % Diffuse Kalman filter
if kalman_algo ~= 4
kalman_algo = 3;
end
[QT,ST] = schur(T);
e1 = abs(ordeig(ST)) > 2-options_.qz_criterium;
[QT,ST] = ordschur(QT,ST,e1);
k = find(abs(ordeig(ST)) > 2-options_.qz_criterium);
nk = length(k);
nk1 = nk+1;
Pinf = zeros(np,np);
Pinf(1:nk,1:nk) = eye(nk);
Pstar = zeros(np,np);
B = QT'*R*Q*R'*QT;
for i=np:-1:nk+2
if ST(i,i-1) == 0
if i == np
c = zeros(np-nk,1);
else
c = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i,i+1:end)')+...
ST(i,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i);
end
q = eye(i-nk)-ST(nk1:i,nk1:i)*ST(i,i);
Pstar(nk1:i,i) = q\(B(nk1:i,i)+c);
Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)';
else
if i == np
c = zeros(np-nk,1);
c1 = zeros(np-nk,1);
else
c = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i,i+1:end)')+...
ST(i,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i)+...
ST(i,i-1)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i-1);
c1 = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i-1,i+1:end)')+...
ST(i-1,i-1)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i-1)+...
ST(i-1,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i);
end
q = [eye(i-nk)-ST(nk1:i,nk1:i)*ST(i,i) -ST(nk1:i,nk1:i)*ST(i,i-1);...
-ST(nk1:i,nk1:i)*ST(i-1,i) eye(i-nk)-ST(nk1:i,nk1:i)*ST(i-1,i-1)];
z = q\[B(nk1:i,i)+c;B(nk1:i,i-1)+c1];
Pstar(nk1:i,i) = z(1:(i-nk));
Pstar(nk1:i,i-1) = z(i-nk+1:end);
Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)';
Pstar(i-1,nk1:i-2) = Pstar(nk1:i-2,i-1)';
i = i - 1;
end
end
if i == nk+2
c = ST(nk+1,:)*(Pstar(:,nk+2:end)*ST(nk1,nk+2:end)')+ST(nk1,nk1)*ST(nk1,nk+2:end)*Pstar(nk+2:end,nk1);
Pstar(nk1,nk1)=(B(nk1,nk1)+c)/(1-ST(nk1,nk1)*ST(nk1,nk1));
end
Z = QT(mf,:);
R1 = QT'*R;
[QQ,RR,EE] = qr(Z*ST(:,1:nk),0);
k = find(abs(diag(RR)) < 1e-8);
if length(k) > 0
k1 = EE(:,k);
dd =ones(nk,1);
dd(k1) = zeros(length(k1),1);
Pinf(1:nk,1:nk) = diag(dd);
end
end
if kalman_algo == 2
no_correlation_flag = 1;
if length(H)==1
H = zeros(nobs,1);
else
if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
H = diag(H);
else
no_correlation_flag = 0;
end