function disc=iterexg(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=rrepmat(Y,[],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);
bind=W;
k=repmat(.5,[1 nrep nval]);
for it=1:nt
	W(it,:,:)=Y(it,:,:).*k.^alph+delta*k;
	C(it,:,:)=exp(g0(it,:,:)+g1(it,:,:).*W(it,:,:)+g2(it,:,:).*W(it,:,:).^2+...
		g3(it,:,:).*W(it,:,:).^3);
	C(it,:,:)=min(C(it,:,:),W(it,:,:));
	k=W(it,:,:)-C(it,:,:);
end
disc=zeros(size(g));
ehat=bet*(C(1:nt-1,:,:)./C(2:nt,:,:)).^gam.*(alph*Y(2:nt,:,:).*(W(1:nt-1,:,:)-C(1:nt-1,:,:)) ...
	.^(alph-1)+delta)-1;
disc(1,:)=sum(sum(ehat));
disc(2,:)=sum(sum(ehat.*W(1:nt-1,:,:)));
disc(3,:)=sum(sum(ehat.*W(1:nt-1,:,:).^2));
if np==4
	disc(4,:)=sum(sum(ehat.*W(1:nt-1,:,:).^3));
end
save iterexd W C ehat