GMM

%=========================================================================
%
% Program to estimate level effect in interest rates by GMM
%
% Code based on Martin, Hurn and Harris, Econometric Time Series Modelling
% Specification, Estimation and Testing
% https://www.cambridge.org/features/econmodelling/chapter10.htm
% 
% This code by Ralf Becker, March 2021
% http://eclr.humanities.manchester.ac.uk/index.php/MATLAB
%=========================================================================


clear all;
clc;
cd 'YOUR DIRECTORY'

% Load data --- monthly December 1946 to February 1991
%     3 month maturity
% extracted from the datafile provided by 
% Martin, Hurn and Harris
% https://www.cambridge.org/features/econmodelling/chapter10.htm

[rt, ~, ~] = xlsread('US3monthRate.xlsx');

drt = trimr(rt,2,0) - trimr(rt,1,1); % creates \Delta r_{t+1}
r1t = trimr(rt,1,1);                 % creates r_t
r2t = trimr(rt,0,2);
t   = length(drt);

%% It is typically good practice to visualise the data
tt = seqa(1946+12/12,1/12,t); % Creates year sequence

subplot(1,2,1);
plot(tt,r1t);
title('r_t')
box off
axis tight

subplot(1,2,2);
plot(tt,drt);
title('\Delta r_{t+1}')
box off
axis tight


%% Estimate the model in the first stage with idendity weighting matrix

ops  = optimset('LargeScale','off','Display','off');
b0   = [0.1;0.1;0.1;1.0];
w0 = eye(length(b0));
bgmm1 = fminunc(@(b) qw(b,drt,r1t,w0),b0,ops);

disp('First Stage GMM estimates')
disp(bgmm1)

%% Now estimate the optimal weighting matrix
% using Newey-West estimator
lmax = 5;   % lag for the NW estimate 
d = meqn(bgmm1,drt,r1t);

% this will calculate Newey-West VCM using lmax lags
s   = d'*d;
tau = 1;
while tau <= lmax
    wtau = d((tau+1):size(d,1),:)'*d(1:(size(d,1)-tau),:);
    s    = s + (1.0-tau/(lmax+1))*(wtau + wtau');
    tau  = tau + 1;
end
w1 = s./t;

% Use this as the weighting matrix for the next pass to the optimisation
% function

%% 2nd Stage 

bgmm2 = fminunc(@(b) qw(b,drt,r1t,w1),bgmm1,ops);

disp('Second Stage GMM estimates')
disp(bgmm2)

%% Further Iterations
% You could run further iterations
% 1) Re-calculate d
% 2) Re-calculate the optimal weighting matrix w based on the new d
% 3) Re-estimate using the new w
%
% For now we stop here
bgmm = bgmm2;
obj = qw(bgmm,drt,r1t,w1);

%% Calculate standard errors
% Compute optimal weigthing matrix at GMM estimates
% using Newey-West estimator
lmax = 5;   % lag for the NW estimate 
d = meqn(bgmm,drt,r1t);

% this will calculate Newey-West VCM using lmax lags
s   = d'*d;
tau = 1;
while tau <= lmax
    wtau = d((tau+1):size(d,1),:)'*d(1:(size(d,1)-tau),:);
    s    = s + (1.0-tau/(lmax+1))*(wtau + wtau');
    tau  = tau + 1;
end
s = s./t;

% Compute standard errors of GMM estimates
dg = numgrad(@meaneqn,bgmm,drt,r1t);
v  = dg'*inv(s)*dg;
cov = inv(v)/t;
se = sqrt(diag(cov));

disp(' ');
disp(['The value of the objective function  = ', num2str(obj) ]);
disp(['J-test                               = ', num2str(t*obj) ]);
disp('Estimates     Std err.   t-stats');
disp( [ bgmm  se  bgmm./se ])
disp(['Newey-West estimator with max lag    = ', num2str(lmax) ]);
disp(' ');


%% Inference t-tests

% Test of gam = 0.0
stat = (bgmm(4) - 0.0)/se(4);
disp(['Test of (gam=0.0) = ', num2str(stat) ]);
disp(['p-value           = ', num2str(2*(1-normcdf(abs(stat)))) ]);
disp(' ');

% Test of gam = 0.5
stat = (bgmm(4) - 0.5)/se(4);
disp(['Test of (gam=0.5) = ', num2str(stat) ]);
disp(['p-value           = ', num2str(2*(1-normcdf(abs(stat)))) ]);
disp(' ');

% Test of gam = 1.0
stat = (bgmm(4) - 1.0)/se(4);
disp(['Test of (gam=1.0) = ', num2str(stat) ]);
disp(['p-value           = ', num2str(2*(1-normcdf(abs(stat)))) ]);
disp(' ');

% Test of gam = 1.5
stat = (bgmm(4) - 1.5)/se(4);
disp(['Test of (gam=1.5) = ', num2str(stat) ]);
disp(['p-value           = ', num2str(2*(1-normcdf(abs(stat)))) ]);
disp(' ');

%% Inference - Overidentifying restrictions

%% Plot volatility function for alternative values of gam
tt = seqa(1946+12/12,1/12,t);
figure(1)

subplot(2,2,1);
plot(tt,drt./r1t.^0.0);
title('$\gamma=0.0$')
box off
axis tight

subplot(2,2,2);
plot(tt,drt./r1t.^0.5);
title('$\gamma=0.5$')
box off
axis tight

subplot(2,2,3);
plot(tt,drt./r1t.^1.0);
title('$\gamma=1.0$')
box off
axis tight

subplot(2,2,4);
plot(tt,drt./r1t.^1.5);
title('$\gamma=1.5$')
box off
axis tight

%
%------------------------- Functions -------------------------------------%
%
%-------------------------------------------------------------------------%
% Define the moment equations 
%-------------------------------------------------------------------------%
function dt = meqn(b,drt,r1t)
    
        ut = drt - b(1) - b(2)*r1t;
        zt = [ones(size(ut,1),1),r1t];
        dt = repmat(ut,1,2).*zt;
        dt = [dt,repmat((ut.^2 - (b(3)^2)*r1t.^(2*b(4)) ),1,2).*zt];
   
end
%-------------------------------------------------------------------------%
% Defines the mean of the moment conditions  
%-------------------------------------------------------------------------%
function ret = meaneqn(b,drt,r1t)

        ret = (mean(meqn(b,drt,r1t)))';

end

%-------------------------------------------------------------------------%
% GMM objective function with  user defined 
% weighting matrix, w
%-------------------------------------------------------------------------%   
function ret = qw(b,drt,r1t,w)
        
    t = length(drt);
    d = meqn(b,drt,r1t);
    g = mean(d)';

    ret = g'*inv(w)*g;
end