% Copyright (C) 2001 Michel Juillard
%
function dr=dr1(iorder,dr,check)

global jacobia_ iy_ ykmin_ ykmax_ gstep_ exo_nbr endo_nbr
global ex_ valf_ it_ exe_ xkmin_ xkmax_ ys_ stdexo_
global fname_ means_ Sigma_e_ lgy_
global eigval options_

xlen = xkmax_ + xkmin_ + 1;
klen = ykmin_ + ykmax_ + 1;
iyv = iy_';
iyv = iyv(:);
iyr0 = find(iyv) ;
it_ = ykmin_ + 1 ;


if exo_nbr == 0
  exe_ = [] ;
end

if ~ iy_(ykmin_+1,:) > 0
  error ('Error in model specification: some variables don"t appear as current') ;
end

if ~check
  if ykmax_ == 0
    error ('No forward variable: no point in using dr1') ;
  end
  if xlen > 1
    error (['SS: stochastic exogenous variables must appear only at the' ...
	    ' current period. Use additional endogenous variables']) ;
  end
end
  
if ykmax_ > 1 & iorder > 1
  error (['Models with leads on more than one period can only be solved' ...
	  ' at order 1'])
end

dr=set_state_space(dr);
kstate = dr.kstate;
kad = dr.kad;
kae = dr.kae;
nstatic = dr.nstatic;
nfwrd = dr.nfwrd;
npred = dr.npred;
nboth = dr.nboth;
order_var = dr.order_var;
nd = size(kstate,1);

sdyn = endo_nbr - nstatic;


tempex = ex_;

it_ = ykmin_ + 1;
z = repmat(dr.ys,1,klen);
z = z(iyr0) ;
%jacobia_=real(diffext('ff1_',[z; exe_])) ;
jacobia_=real(jacob_a('ff1_',[z; exe_])) ;

ex_ = tempex ;
tempex = [];

nz = size(z,1);
k1 = iy_(find([1:klen] ~= ykmin_+1),:);
b = jacobia_(:,iy_(ykmin_+1,order_var));
a = b\jacobia_(:,nonzeros(k1')); 
if any(isinf(a(:)))
  error('DR1: the model doesn''t determine current variables uniquely')
end
if exo_nbr
  fu = b\jacobia_(:,nz+1:end);
end


% buildind D and E
d = zeros(nd,nd) ;
e = d ;

k = find(kstate(:,2) >= ykmin_+2 & kstate(:,3));
d(1:sdyn,k) = a(nstatic+1:end,kstate(k,3)) ;
k1 = find(kstate(:,2) == ykmin_+2);
a1 = eye(sdyn);
e(1:sdyn,k1) =  -a1(:,kstate(k1,1)-nstatic);
k = find(kstate(:,2) <= ykmin_+1 & kstate(:,4));
e(1:sdyn,k) = -a(nstatic+1:end,kstate(k,4)) ;
k2 = find(kstate(:,2) == ykmin_+1);
k2 = k2(~ismember(kstate(k2,1),kstate(k1,1)));
d(1:sdyn,k2) = a1(:,kstate(k2,1)-nstatic);

if ~isempty(kad)
  for j = 1:size(kad,1)
    d(sdyn+j,kad(j)) = 1 ;
    e(sdyn+j,kae(j)) = 1 ;
  end
end

options_ = set_default_option(options_,'qz_criterium',1.000001);
%if ~exist('mjdgges')
if 1
  % using Chris Sim's routines
  [ss,tt,qq,w] = qz(e,d);
  [tt,ss,qq,w] = qzdiv(options_.qz_criterium,tt,ss,qq,w);
  ss1=diag(ss);
  tt1=diag(tt);
  warning_state = warning;
  warning off;
  eigval = ss1./tt1 ;
  warning warning_state;
  nba = nnz(abs(eigval) > options_.qz_criterium);
else
  [ss,tt,w,sdim,eigval,info] = mjdgges(e,d,options_.qz_criterium);
  if info & info ~= nd+2;
    error(['ERROR' info ' in MJDGGES.DLL']);
  end
  nba = nd-sdim;
end

nyf = sum(kstate(:,2) > ykmin_+1);

if check
  dr.rank = rank(w(1:nyf,nd-nyf+1:end));
  dr.eigval = eig(e,d);
  return
end

if nba ~= nyf;
  disp('WARNING: Blanchard-Kahn conditions are not satisfied. Run CHECK to learn more!');
  disp('Press any key to continue');
  pause
end

np = nd - nyf;
n2 = np + 1;
n3 = nyf;
n4 = n3 + 1;
% derivatives with respect to dynamic state variables
% forward variables
gx = -w(1:n3,n2:nd)'\w(n4:nd,n2:nd)';
% predetermined variables
hx = w(1:n3,1:np)'*gx+w(n4:nd,1:np)';
hx = (tt(1:np,1:np)*hx)\(ss(1:np,1:np)*hx);

k1 = find(kstate(n4:nd,2) == ykmin_+1);
k2 = find(kstate(1:n3,2) == ykmin_+2);
dr.ghx = [hx(k1,:); gx(k2(nboth+1:end),:)];

%lead variables actually present in the model
j3 = nonzeros(kstate(:,3));
j4  = find(kstate(:,3));
% derivatives with respect to exogenous variables
if exo_nbr
  a1 = eye(endo_nbr);
  aa1 = [];
  if nstatic > 0
    aa1 = a1(:,1:nstatic);
  end
  dr.ghu = -[aa1 a(:,j3)*gx(j4,1:npred)+a1(:,nstatic+1:nstatic+ ...
						  npred) a1(:,nstatic+npred+1:end)]\fu;
end

% static variables
if nstatic > 0
  temp = -a(1:nstatic,j3)*gx(j4,:)*hx;
  j5 = find(kstate(n4:nd,4));
  temp(:,j5) = temp(:,j5)-a(1:nstatic,nonzeros(kstate(:,4)));
  dr.ghx = [temp; dr.ghx];
  temp = [];
end

if iorder == 1
  return
end

% Second order
tempex = ex_ ;

%hessian = real(hessext('ff1_',[z; exe_]))' ;
kk = flipud(cumsum(flipud(iy_(ykmin_+1:end,order_var)),1));
if ykmin_ > 0
  kk = [cumsum(iy_(1:ykmin_,order_var),1); kk];
end
kk = kk';
kk = find(kk(:));
nk = size(kk,1)+exo_nbr;
k1 = iy_(:,order_var);
k1 = k1';
k1 = k1(:);
k1 = k1(kk);
k2 = find(k1);
kk1(k1(k2)) = k2;
kk1 = [kk1 length(k1)+1:length(k1)+exo_nbr];
kk = reshape([1:nk^2],nk,nk);
kk1 = kk(kk1,kk1);
hessian = zeros(endo_nbr,nk^2);
hessian(:,kk1(:)) = real(jacob2('ff1_',[z; exe_])) ;

ex_ = tempex ;
clear tempex

n1 = 0;
n2 = np;
zx = zeros(np,np);
for i=2:ykmin_+1
  k1 = sum(kstate(:,2) == i);
  zx(n1+1:n1+k1,n2-k1+1:n2)=eye(k1);
  n1 = n1+k1;
  n2 = n2-k1;
end
zx=[zx;dr.ghx;gx*hx;zeros(exo_nbr,np)];
% $$$ k1=nonzeros(iy_(:,order_var)');
% $$$ if exo_nbr 
% $$$   k1=[k1;[length(k1)+1:length(k1)+exo_nbr]'];
% $$$ end
% $$$ kk = reshape([1:length(k1)^2],length(k1),length(k1));  % reorder hessian
% $$$ kk = kk(k1,k1);
% $$$ hessian = hessian(:,kk);
[n1,n2] = size(zx);
if n1*n1*n2*n2 > 1e7
  rhs = zeros(endo_nbr,n2*n2);
  k1 = 1;
  for i1 = 1:n2
      for i2 = 1:n2
	rhs(:,k1) = hessian*kron(zx(:,i1),zx(:,i2));
	k1 = k1 + 1; 
      end
  end
else
  rhs = hessian*kron(zx,zx);
end
rhs = -rhs;

%lhs
hxt = hx';
gx = dr.ghx(end-nyf+1:end,:);
f1 = sparse(jacobia_(:,nonzeros(iy_(ykmin_+2,:))));
% f_y'g_xx hx kron hx + f_y g_xx + f_y'g_x h_xx + f_x' h_xx
%
j1 = sparse(jacobia_(:,iy_(ykmin_+1,order_var)));
% $$$ if ykmin_ > 1;
% $$$   kk = flipud(cumsum(iy_(1:ykmin_,order_var),1))';
% $$$   k1 = find(kk(:,2:ykmin_));
% $$$   k2 = find(kk);
% $$$   kk(k2) = [1:length(k2)];
% $$$   k2 = kk(k1);
% $$$   j1(:,k2) = j1(:,k2)+f1*gx(:,npred+1:end);
% $$$ end
aaa=j1;
j1(:,nstatic+1:nstatic+npred)=...
    j1(:,nstatic+1:nstatic+npred)+f1*gx(:,1:npred);
dr.ghxx = sylvester3(j1,[zeros(endo_nbr,nstatic+npred-nboth) f1],kron(hx,hx),rhs);
%ghxu
%rhs
hu = dr.ghu(nstatic+1:nstatic+npred,:);
zu=[zeros(np,exo_nbr);dr.ghu;gx(:,1:npred)*hu;eye(exo_nbr)];
%kk = reshape([1:np*np],np,np);
%kk = kk(1:npred,1:npred);
%rhs = -hessian*kron(zx,zu)-f1*dr.ghxx(end-nyf+1:end,kk(:))*kron(hx(1:npred,:),hu(1:npred,:));
if n1*n1*n2*exo_nbr > 1e7
  rhs = zeros(endo_nbr,n2*exo_nbr);
  k1 = 1;
  for i1 = 1:n2
      for i2 = 1:exo_nbr
	rhs(:,k1) = hessian*kron(zx(:,i1),zu(:,i2));
	k1 = k1 + 1; 
      end
  end
else
  rhs = hessian*kron(zx,zu);
end
rhs = -rhs-f1*dr.ghxx(end-nyf+1:end,:)*kron(hx,[hu;zeros(np-npred,exo_nbr)]);

%lhs
lhs = j1;
dr.ghxu = lhs\rhs;

%ghuu
%rhs
kk = reshape([1:np*np],np,np);
kk = kk(1:npred,1:npred);
if n1*n1*exo_nbr*exo_nbr > 1e7
  rhs = zeros(endo_nbr,exo_nbr*exo_nbr);
  k1 = 1;
  for i1 = 1:exo_nbr
      for i2 = 1:exo_nbr
	rhs(:,k1) = hessian*kron(zu(:,i1),zu(:,i2));
	k1 = k1 + 1; 
      end
  end
else
  rhs = hessian*kron(zu,zu);
end
rhs = -rhs-f1*dr.ghxx(end-nyf+1:end,kk(:))*kron(hu(1:npred,:),hu(1:npred,:));

%lhs
dr.ghuu = lhs\rhs;


% dr.ghs2
% derivatives of F with respect to forward variables
k1 = reshape([1:nk^2],nk,nk);
k1 = k1(np+endo_nbr+1:np+endo_nbr+nboth+nfwrd,np+endo_nbr+1:np+endo_nbr+nboth+nfwrd);
guu = dr.ghuu(end-nyf+1:end,:);
gu = dr.ghu(end-nyf+1:end,:);
[n1,n2] = size(gu);
if n1*n1*n2*n2 > 1e7
  hessian = hessian(:,k1);
  rhs = zeros(endo_nbr,n2*n2);
  k1 = 1;
  for i1 = 1:n2
      for i2 = 1:n2
	rhs(:,k1) = hessian*kron(gu(:,i1),gu(:,i2));
	k1 = k1 + 1; 
      end
  end
else
  rhs = hessian(:,k1)*kron(gu,gu);
end
rhs = (f1*guu+rhs)*Sigma_e_(:);
dr.fuu = rhs;
rhs = -rhs-dr.fbias;
lhs = j1;
lhs(:,end-nyf+1:end) = lhs(:,end-nyf+1:end) + f1;
dr.ghs2 = lhs\rhs;

% 01/08/2001 MJ put return on iorder == 1 after defining dr.kstate and dr.kdyn
% 01/17/2001 MJ added dr.delta_s: correction factor for order = 2
% 01/21/2001 FC correction of delta_s for more than 1 shock
% 01/23/2001 MJ suppression of redundant sum() in delta_s formula
% 02/22/2001 MJ stderr_ replaced by Sigma_e_
% 08/24/2001 MJ changed the order of the variables, separates static
%               variables and handles only one instance of variables both
%               in lead and lag
% 08/24/2001 MJ added sigma to state variables as in Schmitt-Grohe and
%               Uribe (2001)
% 10/20/2002 MJ corrected lags on several periods bug
% 10/30/2002 MJ corrected lags on several periods bug on static when some
%               intermediary lags are missing
% 12/08/2002 MJ uses sylvester3 to solve for dr.ghxx
% 01/01/2003 MJ added dr.fbias for iterative for dr_algo == 1
% 02/21/2003 MJ corrected bug for models without lagged variables
% 03/02/2003 MJ fixed second order for lag on several periods
% 05/21/2003 MJ add check call argument and make computation for CHECK
% 06/01/2003 MJ added a test for xkmax_ > 1 and order > 1






