This question deals with restricted maximum likelihood (REML) estimation in a particular version of the linear model, namely:
$$ Y = X(\alpha)\beta + \epsilon, \\ \epsilon\sim N_n(0, \Sigma(\alpha)), $$
where $X(\alpha)$ is a ($n \times p$) matrix parametrized by $\alpha \in \mathbb R^k$, as is $\Sigma(\alpha)$. $\beta$ is an unknown vector of nuisance parameters; the interest is in estimating $\alpha$, and we have $k\leq p\ll n$. Estimating the model by maximum likelihood is no problem, but I want to use REML. It is well known, see e.g. LaMotte, that the likelihood $A'Y$, where $A$ is any semi-orthogonal matrix such that $A'X=0$ can be written
$$ L_{\text{REML}}(\alpha\mid Y) \propto\vert X'X\vert^{1/2} \vert \Sigma\vert^{-1/2}\vert X'\Sigma^{-1}X\vert^{-1/2}\exp\left\{-\frac{1}{2} r'\Sigma^{-1}r \right\}, \\ r = (I - X(X'\Sigma^{-1}X)^+X'\Sigma^{-1})Y, $$
when $X$ is full column rank.
My problem is that for some perfectly reasonable, and scientifically interesting, $\alpha$ the matrix $X(\alpha)$ is not of full column rank. All the derivations I have seen of the restricted likelihood above makes use of determinant equalities that are not applicable when $\vert X'X\vert=0$, i.e. they assume full column rank of $X$. This means that the above restricted likelihood is only correct for my setting on parts of the parameter space, and thus is not what I want to optimize.
Question: Are there more general restricted likelihoods, derived, in the statistical literature or elsewhere, without the assumption that $X$ be full column rank? If so, what do they look like?
Some observations:
- Deriving the exponential part is no problem for any $X(\alpha)$ and it may be written in terms of the Moore-Penrose inverse as above
- The columns of $A$ are an (any) orthonormal basis for $C(X)^\bot$
- For known $A$, the likelihood for $A'Y$ can easily be written down for every $\alpha$, but of course the number of basis vectors, i.e. columns, in $A$ depends on the column rank of $X$
If anyone interested in this question believes the exact parameterization of $X,\Sigma$ would help, let me know and I'll write them down. At this point though, I'm mostly interested in a REML for a general $X$ of the correct dimensions.
A more detailed description of the model follows here. Let $y_t = \mu + Ay_{t - 1} + v_t, t = 1, \dots, T$ be an $r$-dimensional first order Vector Autoregression [VAR(1)] where $v_t \overset{iid}{\sim}N(0, \Omega)$. Suppose the process is started in some fixed value $y_0$ at time $t = 0$.
Define $Y = [y_1', \dots, y_T']'$. The model may be written in the linear model form $Y = X\beta + \varepsilon$ using the following definitions and notation:
\begin{align} X &= [1_T \otimes I_r, C^{-1}B] \\ \beta &= [\mu', y_0' - \mu']' \\ \mathrm{var}(\varepsilon)^{-1} &= C'(I_T \otimes \Omega^{-1})C \\ C &= \begin{bmatrix} I_r & 0 & 0 & \cdots \\ -A & I_r & 0 & \cdots \\ 0 & -A & I_r & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{bmatrix} \\ B &= e_{1, T} \otimes A, \end{align}
where $1_T$ denotes a $T-$dimensional vector of ones and $e_{1,T}$ the first standard basis vector of $\mathbb R^T$.
Denote $\alpha = \mathrm{vec}(A)$. Notice that if $A$ is not full rank then $X(\alpha)$ is not full column rank. This includes, for example, cases where one of the components of $y_t$ does not depend on the past.
The idea of estimating VARs using REML is well known in, for example, the predictive regressions literature (see e.g. Phillips and Chen and the references therein.)
It may be worthwhile to clarify that the matrix $X$ is not a design matrix in the usual sense, it just falls out of the model and unless there is a priori knowledge about $A$ there is, as far as I can tell, no way to reparameterize it to be full rank.
I have posted a question on math.stackexchange that is related to this one in the sense that an answer to the math question may help in deriving a likelihood that would answer this question.
