Could the covariance matrix of the moment conditions in GMM be ill-conditioned?

Question

General question:

In a generalized method of moments estimation could the covariance matrix of the moment conditions be ill-conditioned and therefore the inverse not computable?

Background on my model:

I am estimating a random coefficient logit model of demand formulated as a mathematical programm with equilibrium constraints (MPEC) (see Random Coefficients Logit using the MPEC algorithm. There the structural parameters of consumer demand are estimated using simulated GMM. Essentially, I build a model that predicts market shares for products based on the mean valuations of product characteristics (e.g. price, color, size, ...) and the individual deviations from this mean valuation. I then try to find the parameter values that minimize some metric between observed and predicted market shares. Because price is endogenous the metric is not simply the sum of least squares but rather the metric is defined by the instrumental variables and the GMM objective function.

This problem cannot be solved analytically so numerical methods are used. So finding the global minimum requires me to repeatedly solve this problem using different starting values of the parameters I want to estimate. For some values of the starting value I cannot compute the inverse of covariance matrix of moment conditions (efficient estimate of the weighting matrix for two-step GMM). For other starting values this problem does not arise.

Details on my issue:

In my specification for product demand has many indicator variables (300+). Many of these indicators are interactions between indicator variables and are therefore often zero. The matrix of independent variable X (dimension N x K with N observations and K independent variables) also contains "normal" continuous variables but only 5. Therefore the matrix of included and excluded instruments Z (N x L where L is number of instruments) also contains the same amount of zero-one indicator variables.

Now when I try to calculate the covariance-matrix of the moment conditions (to get an efficient estimate of optimal weighting matrix $\hat{W}$) i.e.

$\frac{1}{n} Z'\hat\Omega Z$ with $\hat\Omega = blkdiag(\hat{u}^2_1 ... \hat{u}^2_n )$ and $\hat{u}$ the residuals from the first stage GMM estimate

Matlab will tell me that this matrix is nearly singular and therefore the inverse might be imprecisely calculated.

Warning: Matrix is close to singular or badly scaled. Results may be
inaccurate. RCOND =  3.415469e-18. 

Now note, that Matlab does not tell me hat this matrix is singular. Therefore, the problem is not that I have perfect multicollinaerity due to complete indicators (dummy variable trap).

Is this to be expected with an instrument matrix that includes many indicator variables that are mostly zero for an observation? Therefore is this some kind of numerical problem? Or what else could be the problem?

Answer 1

From what you've told me, I would almost certainly predict that if you: (1) ran standard two stage least squares estimator and (2) tried to compute robust standard errors, then some of the standard errors for your coefficients would either be reported as zero or correctly reported as blank/uncomputable.

My guess at what is happening is that some of your coefficients are exactly identified, hence a heteroskedastic robust covariance matrix cannot be computed for those coefficients. Go to invert that to compute your optimal weighting matrix and Matlab will tell you that the matrix singular or near singular.

Some background

From what you're telling me, Z matrix is full column rank. What you're also telling me is that Z .* U is less than full column rank (where U = u * ones(1, L) and u is your residuals). That is:

$$ \underbrace{Z = \left[\begin{array}{ccc} z_{1,1} & z_{1,2} & z_{z,3} \\ z_{2,1} & z_{2,2} & z_{2,3} \\ \ldots & \ldots & \ldots \end{array} \right]}_{\text{full column rank}} \quad \quad \underbrace{A = \left[\begin{array}{ccc} z_{1,1}u_1 & z_{1,2}u_1 & z_{z,3}u_1 \\ z_{2,1}u_2 & z_{2,2}u_2 & z_{2,3}u_2 \\ \ldots & \ldots & \ldots \end{array} \right]}_{\text{rank deficient?}}$$

The square of the singular values of the matrix A on the right are singular values of $Z'\Omega Z$, that is, $Z' \Omega Z$ has zeros (or near zeros) for singular values if and only if A has zeros or near zeros. (Basically: effective rank = # of singular values above some threshold )

Something I would do is singular value decompositions on both $Z$ and $A$ (or Z'*Z and A'*A, basically same thing for these purposes). Singular values near zero are indicative of a matrix, that from a numerical linear algebra, is near rank deficient. Does $Z$ already have problems without Matlab spitting out an error? Or do the singular values only go to zero (or near zero) once multiplying by the residuals $u_i$?

An educated guess at what might be happening...

Imagine you have 3 observations:

$$ Z = \left[ \begin{array}{cc} 1 & 0 \\ 1 & 1 \\ 1 & 1 \end{array} \right]$$

Z is full column rank because of observation 1. But if the residual $u_1$ for observation 1 is zero...

$$ A = \left[\begin{array}{cc} z_{11}u_1 & z_{12}u_1 \\ z_{21} u_2 & z_{22} u_2 \\ z_{31}u_3 &z_{32}u_3 \end{array}\right] = \left[\begin{array}{cc} 0 & 0 \\ u_2 & u_2 \\ u_3 & u_3 \end{array}\right]$$

Then $A$ is column rank 1 and $Z'\Omega Z = A'A$ will be rank deficient.

Intuitively what would be happening is that if some of your initial b estimates are exactly identified the estimated residual is zero and you can't form an estimate of the variance or covariance! There are 0 leftover degrees of freedom. Mechanically, what gets pumped out is an estimate of zero, then you go to invert that, and Matlab blows up by telling you the matrix is singular or close to singular.

For reference, some MATLAB code I've written/tested to do what I imagine you're trying to do:

Warning, this code does NOT check for rank deficiency. Let endo be your endogenous regressors, exo be your exogenous regressors (constant ones(N, 1) goes here), and instruments be your instruments. Initial setup:

X = [endo, exo];         % N x (L  + K) 
Z = [exo, instruments];  % N x (K + L2)

Calculate your basic ingredients so you don't need to repeat these expensive operations:

ZprimX = Z' * X;
Zprimy = Z' * y;
ZprimZ = Z' * Z;

Get your initial estimates (you will use for weighting matrix):

%%% Method 1 (intuitive): 
% stage1_estimate = Z \ endo;
% endo_hat        = Z * stage1_estimate;   
% b_2sls          = [endo_hat, exo] \ y;
%%% Method 2: (slick)
b_2sls = (Z * (ZprimZ \ ZprimX)) \ y;

Get residuals:

  u_2sls          = y - X * b_2sls;

Calculate weighting matrix inverse

% Basic version
% W_inv = (u_2sls' * u_2sls / N) * ZprimZ / N;

% Robust version
temp = (u_2sls * ones(1, K + L2)) .* Z;
W_inv = temp' * temp / N;

Calculate your estimate b:

%%% Method 1 (intuitive)
% W     = inv(W_inv);
% A     = ZprimX' * W * ZprimX;
% c     = ZprimX' * W * Zprimy;    % = (Zprimy' * W * ZprimX)'
% results.b = inv(A) * c;        

%%% METHOD 2 (slick):
% This is a way without matrix inversion to do it. 
% Solve linear systems rather than invert W_inv
temp = W_inv \ ZprimX;          % Solve inv(W_inv) * ZprimX = W * ZprimX
A    = ZprimX' * temp;
c    = (Zprimy' * temp)';
results.b = A \ c;  

results.b_endo = results.b(1:L);
results.b_exo  = results.b(L+1:end);