function [V]=epstzi9()

% Zeitdiskontfaktor beta < 1/(1+r)
bet=0.96;

% betrachtete Risikoaversionsstufe gamma<>1
gam=2.5;

% Elastizitaet der intertemporalen Substitution als psi<>1
psi=1.2;

% Rechenwerte
egam=1-gam;
emep=1-1/psi;
theta=egam/emep;
emgt=emep;
etheta=emep/egam;
temg=1/emep;

% Anfangsvermoegen
W0=10000;

% Kapitalmarkt und Arbeitseinkommen
load RLdat

% Ueberlebenswahrscheinlichkeiten
load pdat

% maennlich=2, weiblich=3
mw=3;


tm1=cputime;

% letztes Jahr des betrachteten Zeitraums
t=jn;

% fixes Raster fuer die Sparquote s
sW=(0:0.05:0.95);
[sb, sc]=size(sW);

% fixes Raster fuer die Aktienquote a mit gleichverteilter Stuetzstellenzahl
ac=2;
aink=1/(ac-1);
a=(0:aink:1);

% fixes Raster fuer das Vermoegen W
Wink=1000;
W=(Wink:Wink:799*Wink);
[w1c, w2c]=size(W);
Wmax=W(w2c);

% Initialisierung
Vmax=zeros(jn,w2c);
Umax=zeros(jn,w2c);
OptsW=zeros(jn,w2c);
Opta=zeros(jn,w2c);
OptC=zeros(jn,w2c);
Y=zeros(jn,w2c);
sY=zeros(jn,w2c);

SOptsW=zeros(jn,w2c);
SOpta=zeros(jn,w2c);
SOptC=zeros(jn,w2c);
SsY=zeros(jn,w2c);

Vmaxi=zeros(jn,w2c,sn);
Umaxi=zeros(jn,w2c,sn);
OptsWi=zeros(jn,w2c,sn);
Optai=zeros(jn,w2c,sn);
OptCi=zeros(jn,w2c,sn);
Yi=zeros(jn,w2c,sn);
sYi=zeros(jn,w2c,sn);


% Schleife ueber die Anzahl der Simulationen
for sim=1:sn

% Aktienrendite reqt im Jahr t der Simulation
reqt=req(t,sim);

% Arbeitseinkommen Lt im Jahr t der Simulation
Lt=L(t,sim);

% Raster der Konsumwerte C(sW,a,W) im Jahr t
for d1=1:sc
for d2=1:ac
for d3=1:w2c
Cgrid(d1,d2,d3)=(1-sW(d1))*((exp(r+a(d2)*(reqt-r)))*W(d3)+Lt);
end;
end;
end;

Ugrid=((1-bet)*(Cgrid.^emep));

% Ermittlung des aufsummierten Nutzenmaximums Vmax
for d7=1:sc
for d8=1:ac
for d9=1:w2c
Vgrid(d7,d8,d9)=(Ugrid(d7,d8,d9)+bet*1^etheta)^temg;
end;
end;
end;

% Schleife ueber W zur Optimierung von V nach sW und a
for e3=1:w2c

Vsa=[Vgrid(:,:,e3)];

[VsV,IaV]=max(Vsa,[],2);

[V,Is]=max(VsV);

Vmaxi(t,e3,sim)=V;
Umaxi(t,e3,sim)=Ugrid(Is,IaV(Is),e3);
OptsWi(t,e3,sim)=sW(Is);
Optai(t,e3,sim)=a(IaV(Is));
OptCi(t,e3,sim)=Cgrid(Is,IaV(Is),e3);

Yi(t,e3,sim)=(((exp(r+a(IaV(Is))*(reqt-r)))-1)*W(e3)+Lt);
sYi(t,e3,sim)=(Yi(t,e3,sim)-OptCi(t,e3,sim))/Yi(t,e3,sim);

end;

end;

% Potenzieren des Folgenutzens vor der Bildung des Erwartungswerts
EVmaxi=Vmaxi.^egam;

% Schleife ueber W zur Bestimmung des Erwartungswerts von Vmax
for f3=1:w2c

Vmax(t,f3)=mean(EVmaxi(t,f3,:),3);
Umax(t,f3)=mean(Umaxi(t,f3,:),3);
OptsW(t,f3)=mean(OptsWi(t,f3,:),3);
Opta(t,f3)=mean(Optai(t,f3,:),3);
OptC(t,f3)=mean(OptCi(t,f3,:),3);
Y(t,f3)=mean(Yi(t,f3,:),3);
sY(t,f3)=mean(sYi(t,f3,:),3);

SOptsW(t,f3)=std(OptsWi(t,f3,:),0,3);
SOpta(t,f3)=std(Optai(t,f3,:),0,3);
SOptC(t,f3)=std(OptCi(t,f3,:),0,3);
SsY(t,f3)=std(sYi(t,f3,:),0,3);

end;


% naechstes Jahr des betrachteten Zeitraums
t=t-1;


% Schleife ueber alle Vorjahre
while (t>0)

% Schleife ueber die Anzahl der Simulationen
for sim=1:sn

% Aktienrendite reqt im Jahr t der Simulation
reqt=req(t,sim);

% Arbeitseinkommen Lt im Jahr t der Simulation
Lt=L(t,sim);

% Raster der Konsumwerte C(sW,a,W) im Jahr t
for i1=1:sc
for i2=1:ac
for i3=1:w2c
Cgrid(i1,i2,i3)=(1-sW(i1))*((exp(r+a(i2)*(reqt-r)))*W(i3)+Lt);
end;
end;
end;

Ugrid=((1-bet)*(Cgrid.^emep));

% Anlagevermoegen der Folgeperiode
for j1=1:sc
Wt1(j1,:,:)=sW(j1)/(1-sW(j1))*Cgrid(j1,:,:);
end;

% Runden der Wt1-Werte auf das W-Raster zwecks Indexermittlung
% Ermittlung des aufsummierten Nutzenmaximums Vmax
for k1=1:sc
for k2=1:ac
for k3=1:w2c

Wgrid=min(max(Wt1(k1,k2,k3),Wink),Wmax);
Wround=round(Wgrid/Wink)*Wink;
Windex=find(W==Wround);
Vgrid(k1,k2,k3)=(Ugrid(k1,k2,k3)+bet*(p((21+t),mw)*Vmax((t+1),Windex))^etheta)^temg;

end;
end;
end;

% Schleife ueber W zur Optimierung von V nach sW und a
for m3=1:w2c

Vsa=[Vgrid(:,:,m3)];

[VsV,IaV]=max(Vsa,[],2);

[V,Is]=max(VsV);

Vmaxi(t,m3,sim)=V;
Umaxi(t,m3,sim)=Ugrid(Is,IaV(Is),m3);
OptsWi(t,m3,sim)=sW(Is);
Optai(t,m3,sim)=a(IaV(Is));
OptCi(t,m3,sim)=Cgrid(Is,IaV(Is),m3);

Yi(t,m3,sim)=(((exp(r+a(IaV(Is))*(reqt-r)))-1)*W(m3)+Lt);
sYi(t,m3,sim)=(Yi(t,m3)-OptCi(t,m3))/Yi(t,m3);

end;

end;

EVmaxi=Vmaxi.^egam;

% Schleife ueber W zur Bestimmung des Erwartungswerts von Vmax
for n3=1:w2c

Vmax(t,n3)=mean(EVmaxi(t,n3,:),3);
Umax(t,n3)=mean(Umaxi(t,n3,:),3);
OptsW(t,n3)=mean(OptsWi(t,n3,:),3);
Opta(t,n3)=mean(Optai(t,n3,:),3);
OptC(t,n3)=mean(OptCi(t,n3,:),3);
Y(t,n3)=mean(Yi(t,n3,:),3);
sY(t,n3)=mean(sYi(t,n3,:),3);

SOptsW(t,n3)=std(OptsWi(t,n3,:),0,3);
SOpta(t,n3)=std(Optai(t,n3,:),0,3);
SOptC(t,n3)=std(OptCi(t,n3,:),0,3);
SsY(t,n3)=std(sYi(t,n3,:),0,3);

end;


% naechstes betrachtetes Jahr
t=t-1;

end;


Wnow=W0;

for p1=1:jn

Wgrid=min(max(Wnow,Wink),Wmax);
Wnowround=round(Wgrid/Wink)*Wink;
Windex=find(W==Wnowround);

sWm(p1)=OptsW(p1,Windex);
am(p1)=Opta(p1,Windex);
Cm(p1)=OptC(p1,Windex);
Wm(p1)=Cm(p1)/(1-sWm(p1));
Um(p1)=Umax(p1,Windex);
Vm(p1)=Vmax(p1,Windex);
Ym(p1)=Y(p1,Windex);
sYm(p1)=sY(p1,Windex);

sWstd(p1)=SOptsW(p1,Windex);
astd(p1)=SOpta(p1,Windex);
Cstd(p1)=SOptC(p1,Windex);
sYstd(p1)=SsY(p1,Windex);


Wnow=sWm(p1)*Wm(p1);

end;


save xEpstdat sWm am Cm Wm Um Vm Ym sYm sWstd astd Cstd sYstd bet gam psi W0


tm2=cputime-tm1

% Funktionsplot
close all
figure(1)
plot((1:jn),sWm, "linewidth", 4, (1:jn),am, "linewidth", 4);
legend('optimale Sparquote sW', 'optimale Aktienquote a', "location", 'south');
axis([0 (jn+1) -0.1 1.1]);
grid("on");
title('Sparquote sW und Aktienquote a', "fontsize", 18);
xlabel('Periode', "fontsize", 18);
ylabel('optimale Werte', "fontsize", 18);
print -mono -deps -FHelvetica:18 xbilde1.eps

PWlow=min(Wm);
PWhigh=max(Wm);
PClow=min(Cm);
PChigh=max(Cm);
low2=min([PWlow;PClow]);
high2=max([PWhigh;PChigh]);
unten2=low2-0.1*abs(high2-low2);
oben2=high2+0.1*abs(high2-low2);

figure(2)
plot((1:jn),Wm, "linewidth", 4, (1:jn),Cm, "linewidth", 4);
legend('optimales Vermoegen W', 'optimaler Konsum C', "location", 'south');
axis([0 (jn+1) unten2 oben2]);
grid("on");
title('optimales Vermoegen W und optimaler Konsum C', "fontsize", 18);
xlabel('Periode', "fontsize", 18);
ylabel('optimale Werte', "fontsize", 18);
print -mono -deps -FHelvetica:18 xbilde2.eps

PUlow=min(Um);
PUhigh=max(Um);
PVlow=min(Vm);
PVhigh=max(Vm);
low3=min([PUlow;PVlow]);
high3=max([PUhigh;PVhigh]);
unten3=low3-0.1*abs(high3-low3);
oben3=high3+0.1*abs(high3-low3);

figure(3)
plot((1:jn),Vm, "linewidth", 4, (1:jn),Um, "linewidth", 4);
legend('Gesamtnutzen V', 'Nutzen U', "location", 'south');
axis([0 (jn+1) unten3 oben3]);
grid("on");
title('Gesamtnutzen V und Nutzen U', "fontsize", 18);
xlabel('Periode', "fontsize", 18);
ylabel('optimale Werte', "fontsize", 18);
print -mono -deps -FHelvetica:18 xbilde3.eps

b41low=min(sWm);
b41high=max(sWm);
b42low=min(sYm);
b42high=max(sYm);
low4=min([b41low;b42low]);
high4=max([b41high;b42high]);
unten4=low4-0.1*abs(low4);
oben4=high4+0.1*abs(high4);

figure(4)
plot((1:jn),sWm, "linewidth", 4, (1:jn),sYm, "linewidth", 4);
legend('optimale Sparquote sW', 'optimale Sparquote sY', "location", 'south');
axis([0 (jn+1) unten4 oben4]);
grid("on");
title('optimales sW und sY je Periode', "fontsize", 18);
xlabel('Periode', "fontsize", 18);
ylabel('optimale Werte', "fontsize", 18);
print -mono -deps -FHelvetica:18 xbilde4.eps

b51low=min(sWstd);
b51high=max(sWstd);
b52low=min(sYstd);
b52high=max(sYstd);
low5=min([b51low;b52low]);
high5=max([b51high;b52high]);
unten5=low5-0.1*abs(low5);
oben5=high5+0.1*abs(high5);

figure(5)
plot((1:jn),sWstd, "linewidth", 4, (1:jn),sYstd, "linewidth", 4);
legend('Standardabweichung sW', 'Standardabweichung sY', "location", 'south');
axis([0 (jn+1) unten5 oben5]);
grid("on");
title('Standardabweichung je Periode', "fontsize", 18);
xlabel('Periode', "fontsize", 18);
ylabel('std', "fontsize", 18);
print -mono -deps -FHelvetica:18 xbilde5.eps

end;
