%% Expectations
%
%% Notes
% Expectations play a central role in financial economics. For some functions and some random variables, expectations are available in closed forms. However, for most random variables, expectations are not analytically available and so numerical solutions are needed
%
%% Setup
% Clear and reset the workspace and load required data

% Clean up everything
clear all
close all
clc
% Reset rng to make runs the same
rng('default')

%
%% Compute expected values, variances and kurtoses
% Monte carlo uses law of large numbers and i.i.d. shocks to estimate expected values. These examples show that the Monte Carlo values are close to the true value.

% Set number of samples
B = 1000;
% Standard Normal
x = randn(B,1);
ex = mean(x);
ex2 = mean(x.^2);
vx = ex2 - ex.^2;
ex4c = mean((x-ex).^4);
kx = ex4c./vx^2;
disp('Expected Values')
disp([0 1 1 3])
disp('Monte Carlo Values')
disp([ex ex2 vx kx])

% Normal with mean 8% and volatility 20%
x = 0.08 + 0.2 * randn(B,1);
ex = mean(x);
ex2 = mean(x.^2);
vx = ex2 - ex.^2;
ex4c = mean((x-ex).^4);
kx = ex4c./vx^2;
disp('Expected Values')
disp([0.08 .2^2+0.08^2 0.2^2 3])
disp('Monte Carlo Values')
disp([ex ex2 vx kx])

% Students t
v = 8;
x = trnd(v,B,1);
ex = mean(x);
ex2 = mean(x.^2);
vx = ex2 - ex.^2;
ex4c = mean((x-ex).^4);
kx = ex4c./vx^2;
disp('Expected Values')
disp([0 v/(v-2) v/(v-2) (3*v-6)/(v-4)])
disp('Monte Carlo Values')
disp([ex ex2 vx kx])

% Chi2
v = 5;
x = chi2rnd(v,B,1);
ex = mean(x);
ex2 = mean(x.^2);
vx = ex2 - ex.^2;
ex4c = mean((x-ex).^4);
kx = ex4c./vx^2;
disp('Expected Values')
disp([v v^2 + 2*v 2*v 3 + 12/v])
disp('Monte Carlo Values')
disp([ex ex2 vx kx])

%
%% Log-normal
% The expectation of functions other than moments can also be computed using Monte Carlo integration. This example uses the exp of a normal which is, of course, a log-normal.

mu = 0.08;
sigma2 = .2^2;
x = mu + sqrt(sigma2) * randn(B,1);
y = exp(x);

disp('Expected Values')
disp(exp(mu+sigma2/2))
disp('Monte Carlo Values')
disp(mean(y))

%
%% Effect of additional replications
% This example studies the effect of additional draws when using Monte Carlo integration. The second set uses 4 times as many, which produces estimates with a standard deviation 1/2 that of the smaller Monte Carlo (why?)

disp('Expected Values')
disp([v v^2 + 2*v 2*v 3 + 12/v])

mc_values = zeros(10000,4);
for i=1:10000
    v = 5;
    x = chi2rnd(v,B,1);
    ex = mean(x);
    ex2 = mean(x.^2);
    vx = ex2 - ex.^2;
    ex4c = mean((x-ex).^4);
    kx = ex4c./vx^2;
    mc_values(i,:) = ([ex ex2 vx kx]);
end
disp('Std Deviation')
disp(std(mc_values))

mc_values4 = zeros(10000,4);
B4 = B * 4;
for i=1:10000
    v = 5;
    x = chi2rnd(v,B4,1);
    ex = mean(x);
    ex2 = mean(x.^2);
    vx = ex2 - ex.^2;
    ex4c = mean((x-ex).^4);
    kx = ex4c./vx^2;
    mc_values4(i,:) = ([ex ex2 vx kx]);
end
disp('Std Deviation (4x)')
disp(std(mc_values4))
disp('Ratio')
disp(std(mc_values)./std(mc_values4))

%
%% Quadrature
% Quadrature directly attempts to integrate an unknown function by approximating the function with shapes with easy to compute areas. To use quadrature it is necessary to write a function that returns g(x)f(x) where g(x) is the function to take the expectation of (i.e. E[g(x)]) and f(x) is the pdf. In this example, lognormal_quad_target is this function.

mu = 0.08;
sigma = 0.2;
integral_value = integral(@(x)lognormal_quad_target(x,mu,sigma),0.08 - 10 * .2, 0.08 + 10 * .2);
% Older method
quad_value = quadl(@lognormal_quad_target,0.08 - 10 * .2, 0.08 + 10 * .2,[],[],mu,sigma);
disp('Quadrature value')
disp(quad_value)
disp('Integral value')
disp(integral_value)

%
%% Portfolio weights and risk aversion
% This example uses a Taylor expansion of a CRRA utility to compute expected utility which can then be used to find the optimal portfolio allocation. The first section demonstrates one method to compute the required quantities.
% EU\approx\phi_0+\frac{1}{2}\phi_2 E[(r-\mu)^2]+\frac{1}{6}\phi_3 E[(r-\mu)^3] + \frac{1}{24}\phi_4 E[(r-\mu)^4]
% where
% \phi_i = \frac{\partial^i U(W)}{\partial W^i}

W = 1;
gamma = 2;
if gamma ~= 1
    phi_0 = W^(1 - gamma)/(1- gamma);
else
    phi_0 = log(W);
end
phi_1 = W^(-gamma);
phi_2 = -gamma * phi_1 * W^(-1);
phi_3 = (-gamma-1) * phi_2 * W^(-1);
phi_4 = (-gamma-2) * phi_3 * W^(-1);
scales = factorial(0:4).^(-1);

phi = scales .* [phi_0 phi_1 phi_2 phi_3 phi_4];

%
%% Portfolios
% The expected utility can be used to compute the optimal portfolio weights.

nu = 8;
e = trnd(nu,100000,1);
v = nu/(nu-2);
e = e/sqrt(v);
e = 0.08 + 0.2*e;
weights = 0:.02:1;
gammas=1.5:.25:10;
W_0 = 1;

EU = zeros(length(weights),length(gammas));
for j = 1:length(gammas)
    gamma = gammas(j);
    for i = 1:length(weights)
        w = weights(i);
        r = w * e + (1-w) * .01;
        mu = w * 0.08 + (1-w) * 0.01;
        EW_1 = W_0 * (1 + mu);
        mu_prime = mean([ones(size(r)) zeros(size(r)) (r-mu).^2 (r-mu).^3 (r-mu).^4])';
        EU(i,j) = expected_utility(mu_prime,EW_1,gamma);
    end
end

[~,loc] = max(EU);
h=plot(gammas,weights(loc));
set(h,'LineWidth',2)
xlabel('Gamma')
ylabel('Weight on Stock')