Rb at 08:41, 8 February 2022

2022-02-08T08:41:24Z

Rb at 08:40, 8 February 2022

2022-02-08T08:40:35Z

Rb: Created page with "This code illustrates an overidentified (number of moments > number of parameters) 2-step GMM estimation. %=========================================================..."

2022-02-08T08:40:00Z

Created page with "This code illustrates an overidentified (number of moments > number of parameters) 2-step GMM estimation. <source> %=========================================================..."

New page

This code illustrates an overidentified (number of moments > number of parameters) 2-step GMM estimation.

<source>
%=========================================================================
%
% 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
%=========================================================================

clear all;
clc;
cd 'C:\Users\msassrb2\Dropbox (The University of Manchester)\ECON80021\201920\GMM'

% 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(6);
bgmm1 = fminunc(@(b) qw(b,drt,r1t,r2t,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,r2t);

% 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,r2t,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,r2t,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,r2t);

% 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,r2t);
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
J = t*obj;
Jdof = size(s,1)-size(b0,1);
disp(['J-Test of overidentifying restrictions = ', num2str(J) ]);
disp(['p-value = ', num2str(1-chi2cdf(J,2)) ]);
disp(' ');

%% 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,r2t)

ut = drt - b(1) - b(2)*r1t;
zt = [ones(size(ut,1),1),r1t,r2t];
dt = repmat(ut,1,3).*zt;
dt = [dt,repmat((ut.^2 - (b(3)^2)*r1t.^(2*b(4)) ),1,3).*zt];

end
%-------------------------------------------------------------------------%
% Defines the mean of the moment conditions
%-------------------------------------------------------------------------%
function ret = meaneqn(b,drt,r1t,r2t)

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

end

%-------------------------------------------------------------------------%
% GMM objective function with user defined
% weighting matrix, w
%-------------------------------------------------------------------------%
function ret = qw(b,drt,r1t,r2t,w)

t = length(drt);
d = meqn(b,drt,r1t,r2t);
g = mean(d)';

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

</source>

@@ Line 1: / Line 1: @@
 This code illustrates an overidentified (number of moments > number of parameters) 2-step GMM estimation.
+This code requires numgrad.m in your working directory. This is a function to calculate numerical gradients.
 <source>

@@ Line 15: / Line 15: @@
 clear all;
 clc;
-cd 'C:\Users\msassrb2\Dropbox (The University of Manchester)\ECON80021\201920\GMM'
+cd 'C:\Users\YOUR_WORKING_DIRECTORY'
 % Load data --- monthly December 1946 to February 1991

GMM over - Revision history

Rb at 08:41, 8 February 2022

Rb at 08:40, 8 February 2022

Rb: Created page with "This code illustrates an overidentified (number of moments > number of parameters) 2-step GMM estimation. %=========================================================..."