%% Method of Moments
%
%% Notes
% This solution relies heavily on the maximum likelihood solution, especially for the Probit standard error Monte Carlo.
% This solution first shows a simple example, and then wraps it in a loop to perform the Monte Carlo.
%
%% 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')
% Load the data
load SP_FTSE

%
%% S&P 500
% This solution answers questions 1 - 4 for the S&P 500 and then for the FTSE 100. This produces a simpler solution.
%
%% Returns
% Here I use 100 x returns which are nicer to look at than actual returns. The method of moment estiamtors are implemented directly from the notes.

% 100 x returns
r = 100 * diff(log(SP500_monthly(:,3)));
% First four moments, and transformation to get skewness and kurtotis
mu = mean(r)
sigma2 = mean((r-mu).^2)
mom3 = mean((r-mu).^3);
skewness = mom3/sigma2^(3/2)
mom4 = mean((r-mu).^4);
kurtosis = mom4/sigma2^2

%
%% Covariance estimation
% Here I use compute the moment errors for the mean and the variance. Since these conditions are linear in the moments I am interested in, the Jacobian is just an identity matrix and so it doesn't appear.

u1 = r - mu;
u2 = (r - mu).^2 - sigma2;
u = [u1 u2];
T = length(r);
% This is the cov that appears in the sqrt(T)(theta-hat - theta_0) version
% of the CLT
acov = u'*u / T
% This is the covariance uses for inference, ntoe the divide by T
cov_for_inverence = acov / T
% The standard error is just the square root of the parameter variances
stderr = sqrt(diag(cov_for_inverence))

%
%% Sharpe ratio
% The sharpe ratio is the ratio of the mean to the standard deviation. The delta method is used to compute the covariance of the Sharpe ratio

SR = mu/sqrt(sigma2)
% Derivative matrix of the non-linear function w.r.t. parameters
D = [sigma2^(-1/2) -1/2*mu*sigma2^(-3/2)];
% Asymptotic covariance
SR_acov = D*acov*D'
% Divide by T for inference
SR_cov_for_inference = SR_acov/T
% Standard error
SR_stderr = sqrt(diag(SR_cov_for_inference))
SR_tstat = SR/SR_stderr

%
%% Moving windows
% Moving windows are simple graphical tools to examine estimators for stability. This code uses the same moment estimators are in the first part of the solution. subplot is used to graph the data.

% Relabel returns so I can use r in the loop
returns = r;
% Initialize somewhere to store the estimates
mus = zeros(T-119,1);
sigma2s = zeros(T-119,1);
skewnesses = zeros(T-119,1);
kurtoses = zeros(T-119,1);
% Use a loop to estimate the parameters
for i=1:T-119
    r = returns(i:i+119);
    mus(i) = mean(r);
    sigma2s(i) = mean((r-mus(i)).^2);
    mom3 = mean((r-mus(i)).^3);
    skewnesses(i) = mom3/sigma2s(i)^(3/2);
    mom4 = mean((r-mus(i)).^4);
    kurtoses(i) = mom4/sigma2s(i)^2;
end
% Form into a matrix
data = [mus sigma2s skewnesses kurtoses];
% Open a new figure
figure(1)
% Loop over data and use subplot to graph
for i=1:4
    subplot(4,1,i)
    plot(data(:,i))
    axis tight
end

%
%% FTSE 100
% The remainder presents the results for the FTSE. Note the differences in teh kurtosis, Sharpe ratio and the standard errors
%
%% Returns
% Here I use 100 x returns which are nicer to look at than actual returns. The method of moment estiamtors are implemented directly from the notes.

% 100 x returns
r = 100 * diff(log(FTSE100_monthly(:,3)));

% First four moments, and transformation to get skewness and kurtotis
mu = mean(r)
sigma2 = mean((r-mu).^2)
mom3 = mean((r-mu).^3);
skewness = mom3/sigma2^(3/2)
mom4 = mean((r-mu).^4);
kurtosis = mom4/sigma2^2

%
%% Covariance estimation
% Here I use compute the moment errors for the mean and the variance. Since these conditions are linear in the moments I am interested in, the Jacobian is just an identity matrix and so it doesn't appear.

u1 = r - mu;
u2 = (r - mu).^2 - sigma2;
u = [u1 u2];
T = length(r);
% This is the cov that appears in the sqrt(T)(theta-hat - theta_0) version
% of the CLT
acov = u'*u / T
% This is the covariance uses for inference, ntoe the divide by T
cov_for_inverence = acov / T
% The standard error is just the square root of the parameter variances
stderr = sqrt(diag(cov_for_inverence))

%
%% Sharpe ratio
% The sharpe ratio is the ratio of the mean to the standard deviation. The delta method is used to compute the covariance of the Sharpe ratio

SR = mu/sqrt(sigma2)
% Derivative matrix of the non-linear function w.r.t. parameters
D = [sigma2^(-1/2) -1/2*mu*sigma2^(-3/2)];
% Asymptotic covariance
SR_acov = D*acov*D'
% Divide by T for inference
SR_cov_for_inference = SR_acov/T
% Standard error
SR_stderr = sqrt(diag(SR_cov_for_inference))
SR_tstat = SR/SR_stderr

%
%% Moving windows
% Moving windows are simple graphical tools to examine estimators for stability. This code uses the same moment estimators are in the first part of the solution. subplot is used to graph the data.

% Relabel returns so I can use r in the loop
returns = r;
% Initialize somewhere to store the estimates
mus = zeros(T-119,1);
sigma2s = zeros(T-119,1);
skewnesses = zeros(T-119,1);
kurtoses = zeros(T-119,1);
% Use a loop to estimate the parameters
for i=1:T-119
    r = returns(i:i+119);
    mus(i) = mean(r);
    sigma2s(i) = mean((r-mus(i)).^2);
    mom3 = mean((r-mus(i)).^3);
    skewnesses(i) = mom3/sigma2s(i)^(3/2);
    mom4 = mean((r-mus(i)).^4);
    kurtoses(i) = mom4/sigma2s(i)^2;
end
% Form into a matrix
titles = {'Ann mean','Ann std','Skewness','Kurtosis'};
data = [12 * mus / 100 sqrt(12*sigma2s/10000) skewnesses kurtoses];
% Open a new figure
figure(2)
% Loop over data and use subplot to graph
for i=1:4
    subplot(4,1,i)
    plot(data(:,i))
    title(titles{i});
    axis tight
end