function [gx,hx] = gxhx(fy,fx,fyp,fxp,stake);
%This program computes the matrices gx and hx that define the first-order approximation 
%of the DSGE model. That is, if 
%E_t[f(yp,y,xp,x)=0, then the solution is of the form
%xp = h(x,sigma) + sigma * eta * ep
%y = g(x,sigma).
%The first-order approximations to the functions g and h around the point (x,sigma)=(xbar,0), where xbar=h(xbar,0), are:
%h(x,sigma) = xbar + hx (x-xbar) 
%and
%g(x,sigma) = ybar + gx * (x-xbar),
%where ybar=g(xbar,0). 
%Inputs: fy fyp fx fxp
%Outputs: gx hx
%(c) Stephanie Schmitt-Grohe and Martin Uribe
%Date: February 11, 2005

if nargin<5
    stake = 1;
end


A = [-fxp -fyp];
B = [fx fy];

[V,D]=eig(B,A,'qz');
cs=find(abs(diag(D))<stake); %Find the columns corresponding to the roots of D with modulus < stake
ncs=numel(cs); 
NX = size(fx,2);
if ncs<NX
    warning('no local equilirium')
elseif ncs>NX
    warning('eqm is indeterminate')
end

la = D(cs,cs);
Vs = V(:, cs);
P = Vs(1:ncs,:);

if rank(P)<ncs;
	error('Invertibility condition violated')
end

hx =  P*la/P;
gx = Vs(ncs+1:end,:)/P;
hx=real(hx);
gx=real(gx);