%% FOCS for the stochastic growth model

function equ = focsSGM(park, Grid,fspace);
 
global alpha betta sig rho delta A_bar
% global LowerBound UpperBound
global nQuadr QuadrWeights QuadrPoints

LowerBound = fspace.a;
UpperBound = fspace.b;


%rename grid
k = Grid(:,1);
A = Grid(:,2);

% evaluate policy functions
knext = funeval(park,fspace, Grid);
fofk = exp(A).*(k.^alpha) + (1-delta).*k;
% c = fofk - knext;
c = max(fofk - knext, zeros(length(Grid),1));


n = length(k);
% generate nQuadr replications of the Grid, one for each realization of shock:
Grid_knext = kron(knext,ones(nQuadr,1));

% Expected value of next period A, corresponding to Grid:
ExpA =  rho*A;
% all realizations of next A:
GridANext = kron(ExpA,ones(nQuadr,1)) + ...
    kron(ones(n,1),QuadrPoints);
% truncate it to state space:
GridANext = min(max(GridANext,LowerBound(2)),UpperBound(2));

GridNext = [Grid_knext GridANext];


% calculate variables at t+1
knextnext = funeval(park,fspace, GridNext);
fofknext = exp(GridANext).*(Grid_knext.^alpha) + (1-delta).*Grid_knext;
 cnext = max(fofknext - knextnext, zeros(length(Grid_knext),1));
% cnext = fofknext - knextnext;
mucnext = muc(cnext);
mpknext = mpk(GridNext);
% calculate expectations with quadrature
exp_mucnext = (QuadrWeights'*reshape(mpknext.*mucnext,nQuadr,n))';





% equation to be solved: Euler equation
equ = (muc(c) - betta.*exp_mucnext)./muc(c);

% avoid strange solutions
if (any(cnext<0))  || (any(c<0)) || (any(knext<0))  || (any(knextnext<0)) 
    equ(1) = 1e100;
end;

