%GXX_HXX_SPARSE.M
%[gxx,hxx] = gxx_hxx_sparse(fx,fxp,fy,fyp,fypyp,fypy,fypxp,fypx,fyyp,fyy,fyxp,fyx,fxpyp,fxpy,fxpxp,fxpx,fxyp,fxy,fxxp,fxx,hx,gx) 
%finds the 3-dimensional arrays gxx and hxx necessary to compute the 2nd order approximation 
%to the decision rules of a DSGE model of the form E_tf(yp,y,xp,x)=0, with solution 
%xp=h(x,sigma) + sigma * eta * ep and y=g(x,sigma). For more details, see  
%``Solving Dynamic General Equilibrium Models Using a Second-Order Approximation to the Policy 
%Function,'' by Stephanie Schmitt-Grohe and Martin Uribe, JEDC, January 2004, p. 755-775. 
%
%INPUTS: First and second derivatives of f and first-order approximation to the functions g and 
%h: fx, fxp, fy, fyp, fypyp, fypy, fypxp, fypx, fyyp, fyy, fyxp, fyx, fxpyp, fxpy, fxpxp, fxpx, 
%fxyp, fxy, fxxp, fxx, hx, gx
%
%OUTPUTS: Second-order derivatives of the functions g and h with respect to x, evaluated 
%at (x,sigma)=(xbar,0), where xbar=h(xbar,0). That is, hxx gxx
%
% We solve a linear system of the type q = Q * x where x is a vector containing the elements of 
%gxx and hxx appropritely stacked and q and Q are, respectively, a vector and a matrix whose 
%elements are functions of the inputs of the program. 
%
%(c) Stephanie Schmitt-Grohe and Martin Uribe
%
%Date February 18, 2004, modified 2/20/2005 to make matrices sparse.
%Difference with gxx_hxx.m:  matrices are sparse. This allows one to handle larger economies.
%gxx_hxx_sparse.m is slower than gxx_hxx.m, but runs less often into memory problems 

function [gxx,hxx] = gxx_hxx_sparse(fx,fxp,fy,fyp,fypyp,fypy,fypxp,fypx,fyyp,fyy,fyxp,fyx,fxpyp,fxpy,fxpxp,fxpx,fxyp,fxy,fxxp,fxx,hx,gx)

m=0;
nx = size(hx,1); %rows of hx and hxx
ny = size(gx,1); %rows of gx and gxx
n = nx + ny; %length of f
ngxx = nx^2*ny; %elements of gxx

sg = [ny nx nx]; %size of gxx
sh = [nx nx nx]; %size of hxx

%Q = zeros(n*nx*(nx+1)/2,n*nx*nx);
%q = zeros(n*nx*(nx+1)/2,1);
Q = sparse(n*nx*(nx+1)/2,n*nx*nx);
q = sparse(n*nx*(nx+1)/2,1);

gxx=zeros(sg);
hxx=zeros(sh);
GXX=zeros(sg);
HXX=zeros(sh);


for i=1:n
for j=1:nx
%for k=1:nx
for k=1:j
m = m+1;

%First Term
q(m,1) = ( shiftdim(fypyp(i,:,:),1) * gx * hx(:,k) + shiftdim(fypy(i,:,:),1) * gx(:,k) + shiftdim(fypxp(i,:,:),1) *  hx(:,k) + shiftdim(fypx(i,:,k),1) )' * gx * hx(:,j); 

% Second term

GXX(:) = kron(ones(nx^2,1),fyp(i,:)');

pGXX = permute(GXX,[2 3 1]);
pGXX(:) = pGXX(:) .* kron(ones(nx*ny,1),hx(:,j));
GXX=ipermute(pGXX,[2 3 1]);

pGXX = permute(GXX,[3 1 2]);
pGXX(:) = pGXX(:) .* kron(ones(nx*ny,1),hx(:,k));
GXX=ipermute(pGXX,[3 1 2]);

Q(m,1:ngxx)=GXX(:)';

GXX=0*GXX;

%Third term

HXX(:,j,k) = (fyp(i,:) * gx)';

Q(m,ngxx+1:end)=HXX(:)';

HXX = 0*HXX;

%Fourth Term
q(m,1) = q(m,1) + ( shiftdim(fyyp(i,:,:),1) * gx * hx(:,k) +  shiftdim(fyy(i,:,:),1) * gx(:,k) + shiftdim(fyxp(i,:,:),1) * hx(:,k) +  shiftdim(fyx(i,:,k),1) )' * gx(:,j); 

%Fifth Term

GXX(:,j,k)=fy(i,:)';

Q(m,1:ngxx) = Q(m,1:ngxx) + GXX(:)';

GXX = 0*GXX;

%Sixth term
q(m,1) = q(m,1) + ( shiftdim(fxpyp(i,:,:),1) * gx * hx(:,k) + shiftdim(fxpy(i,:,:),1) * gx(:,k) + shiftdim(fxpxp(i,:,:),1) * hx(:,k) + fxpx(i,:,k)')' * hx(:,j);


%Seventh Term

HXX(:,j,k)=fxp(i,:)';

Q(m,ngxx+1:end) = Q(m,ngxx+1:end) + HXX(:)';

HXX = 0*HXX;

%Eighth Term
q(m,1) = q(m,1) +  shiftdim(fxyp(i,j,:),1) * gx * hx(:,k) +  shiftdim(fxy(i,j,:),1) * gx(:,k) +  shiftdim(fxxp(i,j,:),1) * hx(:,k) + fxx(i,j,k);

end %k 
end %j
end %i

A = temp(nx,ny);

Qt = Q* A;

xt = -Qt\q;
x = A* xt;
x=full(x);
gxx(:)=x(1:ngxx);
hxx(:) = x(ngxx+1:end);


function A = temp(nx,ny) %subfunction
%function A = temp(nx,ny)
%This function creates a matrix A, size n*nx^2 by n*nx*(nx+1)/2, such that x = A xtilde, 
%where x is the vector  [gxx(:); hxx(:)] and xtilde is a vector containing
%all the elements  gxx(i,j,k) and hxx(i,j,k), respectively, such that j<=
%k, that is, it is a subset of the elements in the vector x that appear only once.
%This is so because gxx(i,j,k) is symmetric with respect to j and k, that is, gxx(i,j,k)=gxx(i,k,j).
%The reason we use this program is that the matrix we invert (which in the earlier versioni of this program used to be
%Q and now is Qt) is of  size n*nx*(nx+1)/2 by n*nx*(nx+1)/2 rather than of size
%n*nx^2 by n*nx^2. This reduces computation time by about 30 percent. (Feb 16, 2004)
%(April 20, 2005) Made A a sparse matrix. This allows us to handle cases with larger values of ny and nx. 

Ahxx=sparse(nx^3, nx^2*(nx+1)/2);
Agxx=sparse(ny*nx^2, ny*nx*(nx+1)/2);
mx=0;
my=0;
for k=1:nx;
    for j=k:nx;
        for i=1:nx; 
            mx=mx+1;
            Ahxx((j-1)*nx+i+(k-1)*nx*nx, mx) = 1;
            Ahxx((k-1)*nx+i+(j-1)*nx*nx, mx) = 1;
        end %i=1:nx
        for i=1:ny; 
            my=my+1;
            Agxx((j-1)*ny+i+(k-1)*ny*nx, my) = 1;
            Agxx((k-1)*ny+i+(j-1)*ny*nx, my) = 1;
        end %i=1:ny
    end %j
end %k
[ig,jg,sg] = find(Agxx);
[ih,jh,sh]=find(Ahxx);
ih=ih+ny*nx^2;
jh=jh+ny*nx*(nx+1)/2;
A = sparse([ig; ih], [jg; jh], [sg; sh],(nx+ny)*nx^2, (nx+ny)*nx*(nx+1)/2); 
%A(1:ny*nx^2, 1:ny*nx*(nx+1)/2)=Agxx;
%A(ny*nx^2+1:end, ny*nx*(nx+1)/2+1:end)=Ahxx;