function disc=glerkg(g,Y)
%function disc=iterexg(g,Y)
bet=.95; alph=.3; gam=2; delta=1-.07;
[np,nval]=size(g);
[nt,nrep]=size(Y);
Y=repmat(Y,[1 1 nval]);
g0=rrepmat(g(1,:)',[nt,nrep],1);
g1=rrepmat(g(2,:)',[nt,nrep],1);
g2=rrepmat(g(3,:)',[nt,nrep],1);
if np==4
	g3=rrepmat(g(4,:)',[nt,nrep],1);
else
	g3=zeros(size(g2));
end
C=zeros(nt,nrep,nval);
%W=zeros(nt,nrep,nval);
W=.5:.1:7;
W=repmat(W',[1 nrep nval]);
C=exp(g0+g1.*W+g2.*W.^2+g3.*W.^3);
C=min(C,W);
k=W-C;
Wnew=Y.*k.^alph+delta*k;
Cnew=exp(g0+g1.*Wnew+g2.*Wnew.^2+g3.*Wnew.^3);
Cnew=min(Cnew,Wnew);
disc=zeros(size(g));
ehat=zeros(size(W));
ndx=find(C>0&k>0);
ehat(ndx)=bet*(Cnew(ndx)./C(ndx)).^(-gam).*(alph*Y(ndx).*k(ndx).^(alph-1)+delta)-1;
ndx=find(C==0|k==0);
ehat(ndx)=1e10;
disc(1,:)=sum(sum(ehat));
disc(2,:)=sum(sum(ehat.*(W-3.5)));
disc(3,:)=sum(sum(ehat.*(W.^2+8.5-W*7.1)));
if np==4
	disc(4,:)=sum(sum(ehat.*W.^3));
end
save glerkg W C ehat g