% Macro 3. Econ202C. Spring 2006. Homework 1. Problem 2. Anton Cheremukhin
clear all;
close all;
clc;
% get data
[time, re, dlc, rf] = textread('EPData.txt','%d %f %f %f');

% calculate moments
eret=(re'*ones(size(re)))/size(re,1)
ester=sqrt(((re-eret)'*(re-eret))/size(re,1))
rret=(rf'*ones(size(re)))/size(rf,1)
rster=sqrt(((rf-rret)'*(rf-rret))/size(rf,1))
cret=(dlc'*ones(size(re)))/size(dlc,1);
cster=sqrt(((dlc-cret)'*(dlc-cret))/size(dlc,1));

ccorr=(((dlc(1:(size(dlc,1)-1),1)-cret)'*(dlc(2:size(dlc,1),1)-cret))/((size(dlc,1)-1)*(cster^2)));

%calibrate parameters
mu=cret
delta=cster
fi=(ccorr+1)/2;
lambda=[1+mu-delta 1+mu+delta]
epremium=eret-rret

beta0=[0.9 0.925 0.95 0.975 0.98 0.99 0.995 0.999];
alpha0=1:0.5:20;
for i=1:size(beta0,2)
    for j=1:size(alpha0,2)
beta=beta0(i);
alpha=alpha0(j);
%beta=0.995;
%alpha=3;
P=[fi 1-fi
   1-fi fi];
pi=[1/2 1/2];
Rf=pi*(1./(beta*P*(lambda'.^(-alpha)))-1);
%A=eye(2,2)-beta*P.*((lambda'.^(1-alpha))*ones(1,2));
A=eye(2,2)-beta*P*diag(lambda.^(1-alpha));
B=beta*P*(lambda'.^(1-alpha));
w=inv(A)*B;
r=(ones(2,1)*(lambda'.*(w+1))')./(w*ones(1,2))-1;
Re=pi*(P.*r)*ones(2,1);
dr=Re-Rf;
rf0((i-1)*size(alpha0,2)+j)=Rf;
dr0((i-1)*size(alpha0,2)+j)=dr;

    end
end
figure(1)
plot(rf0,dr0,'.',rret,epremium,'or',...
               'LineWidth',2)
xlabel('risk-free rate');
ylabel('equity premium');
legend('simulations','data',2)

%HJ bounds
for i=1:size(beta0,2)
    for j=1:size(alpha0,2)
%        beta=0.995;
%        alpha=1.5;
    beta=beta0(i);
    alpha=alpha0(j);
    y=beta*((1./(1+dlc)).^alpha);
    k=(i-1)*size(alpha0,2)+j;
    Ey(k)=y'*ones(size(y))/size(y,1);
    Vy(k)=((y-Ey(k))'*(y-Ey(k)))/size(y,1);
    end
end

Er=0.75:0.02:0.99;
m=size(Er,2);
Ve=(((1-(1+eret)*Er).^2)/ester).^2;
Vrf=(((1-(1+rret)*Er).^2)/rster).^2;
figure(2)
plot(sqrt(Ve),Er,'-',sqrt(Vrf(m-7:m)),Er(m-7:m),'--',sqrt(Vy),Ey,'.',...
               'LineWidth',2)
xlabel('standard error');
ylabel('expected return');
legend('HJBound for equity','HJBound for bond','simulations')