Difference between revisions of "GMM over"
(Created page with "This code illustrates an overidentified (number of moments > number of parameters) 2-step GMM estimation. <source> %=========================================================...") |
|||
Line 15: | Line 15: | ||
clear all; | clear all; | ||
clc; | clc; | ||
− | cd 'C:\Users\ | + | cd 'C:\Users\YOUR_WORKING_DIRECTORY' |
% Load data --- monthly December 1946 to February 1991 | % Load data --- monthly December 1946 to February 1991 |
Revision as of 08:40, 8 February 2022
This code illustrates an overidentified (number of moments > number of parameters) 2-step GMM estimation.
%=========================================================================
%
% 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\YOUR_WORKING_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(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