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

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?

• I take it you are doing IV from a GMM perspective and in Stata lingo, trying to calculate the "robust" weighting matrix? What does Matlab think is the rank of Z'*Z (i.e. rank(Z'*Z))? What is the smallest singular value in S if you do [U, S, V] = svd(Z'*Z)? – Matthew Gunn May 31 '16 at 21:48
• Yes, I am using GMM to estimate IV with hetereoskedasticity robust standard error. Z'*Z is what I use for the first stage of GMM to calculate a consistent estimate of weighting matrix. The rank of Z'Z is 361 = L and therefore of full rank. But the problem does not arise until the second stage where I estimate $\frac{1}{n}Z'\hat{\Omega}Z$. I'll see that I get the rank there too. – TiredHornet Jun 1 '16 at 12:51
• Are a bunch of the residuals $u_i$ zero or extremely close to zero? (eg like e-14)? And how big is the smallest singular value of $Z'*Z$? (i.e. is it already close to being rank deficient)? – Matthew Gunn Jun 1 '16 at 13:30
• Let $u$ be the residuals from your initial estimate. What is sum(abs(u) < e-13)? If the coefficients on one or more of your indicator variables are exactly identified, you won't be able to make heteroskedastic robust estimate of your initial estimator and hence can't form the weighting matrix. See my answer... – Matthew Gunn Jun 1 '16 at 20:20
• I ran [U, S, V] = svd(Z'*Z) and min(diag(S)) = 0.702. I am still waiting on the residuals. It takes many hours to solve the random coefficient logit model of which the GMM estimation is a part. I'll get back to you when I have them. – TiredHornet Jun 1 '16 at 20:57

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;


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

• Bonus points for showing alternative to explicit matrix inversion. – Mark L. Stone May 31 '16 at 22:16
• @matthew-gunn Thank you for your extensive examples on GMM. The nice solutions without using inv(). I have some similar solutions myself. But I actually need to calculate the weighting matrix explicitly as this GMM problem is part of a random coefficent logit model of demand formulated as an MPEC. I need the weighting matrix to the calculate the hessian of the constraints. Complete indicators should not be a problem. Matlab would otherwise give me a different error and tell me that the matrix I am trying to invert is in fact singular. I have added this not to my question. – TiredHornet Jun 1 '16 at 12:53
• @TiredHornet Is "MATLAB Code for BLP’s GMM estimator of the Random Coefficients Logit using the MPEC algorithm" faculty.chicagobooth.edu/jean-pierre.dube/research/… sort of what you're doing? – Mark L. Stone Jun 1 '16 at 13:18
• @mark-l-stone. Yes my code is based on that with adjustments to my specification and using sparse-grid intergration instead of monte-carlo integration. – TiredHornet Jun 1 '16 at 19:51
• I did as you advised. The smallest singular value of $Z′∗Z$ is 0.7264 while the smallest of $A′∗A$ is 1.4437e-10. So it seems like I have the problem you describe above. I have to read up on SVD since I never came into contact with that before. But is there a way to find out which variable is causing the issue? (had to repost due to mistakes in perceive versions of the comment) – TiredHornet Jun 2 '16 at 14:19