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