Consider the following estimator
\[\hat{\boldsymbol{\theta}} = \operatorname{argmin}_{\boldsymbol{\theta}} \; \left(\mathbf{X}\boldsymbol{\theta} - \hat{\boldsymbol{\nu}}\right)^T \boldsymbol{\Omega} \left(\mathbf{X}\boldsymbol{\theta} - \hat{\boldsymbol{\nu}}\right),\]
subject to
\[\mathbf{X} \hat{\boldsymbol{\theta}} \geq \hat{\boldsymbol{\nu}}.\]
Assuming \(\boldsymbol{\Omega}\) to be a diagonal matrix (with positive entries) we have
\[\left(\mathbf{X}\boldsymbol{\theta} - \hat{\boldsymbol{\nu}}\right)^T \boldsymbol{\Omega} \left(\mathbf{X}\boldsymbol{\theta} - \hat{\boldsymbol{\nu}}\right) = \left(\boldsymbol{\Omega}^{1/2} \mathbf{X}\boldsymbol{\theta} - \boldsymbol{\Omega}^{1/2} \hat{\boldsymbol{\nu}}\right)^T \left(\boldsymbol{\Omega}^{1/2} \mathbf{X}\boldsymbol{\theta} - \boldsymbol{\Omega}^{1/2} \hat{\boldsymbol{\nu}}\right) = \left(\mathbf{X}^*\boldsymbol{\theta} - \hat{\boldsymbol{\nu}}^*\right)^T \left(\mathbf{X}^*\boldsymbol{\theta} - \hat{\boldsymbol{\nu}}^*\right),\] where \(\mathbf{X}^* = \boldsymbol{\Omega}^{1/2} \mathbf{X}\) and \(\hat{\boldsymbol{\nu}}^* = \boldsymbol{\Omega}^{1/2} \hat{\boldsymbol{\nu}}\). Then, using the results of Liew, 1976, we have
where \(\tilde {\boldsymbol{\theta}}\) corresponds to the standard GMWM estimator, i.e.
\[\begin{align} \hat{\boldsymbol{\theta}} &= \tilde {\boldsymbol{\theta}} + \left[\left(\mathbf{X}^*\right)^{T} \mathbf{X}^*\right]^{-1} \mathbf{X}^T \left\{\mathbf{X} \left[\left(\mathbf{X}^*\right)^{T} \mathbf{X}^*\right]^{-1}\mathbf{X}^T\right\}^{-1} \left(\hat{\boldsymbol{\nu}} - \mathbf{X} \hat{\boldsymbol{\theta}}\right)\\ &= \tilde {\boldsymbol{\theta}} + \left(\mathbf{X}^{T} \boldsymbol{\Omega} \mathbf{X}\right)^{-1} \mathbf{X}^T \left\{\mathbf{X} \left(\mathbf{X}^{T} \boldsymbol{\Omega} \mathbf{X}\right)^{-1}\mathbf{X}^T\right\}^{-1} \left(\hat{\boldsymbol{\nu}} - \mathbf{X} \hat{\boldsymbol{\theta}}\right)\\ \end{align}\]
\[\tilde {\boldsymbol{\theta}} = \operatorname{argmin}_{\boldsymbol{\theta}} \; \left(\mathbf{X}\boldsymbol{\theta} - \hat{\boldsymbol{\nu}}\right)^T \boldsymbol{\Omega} \left(\mathbf{X}\boldsymbol{\theta} - \hat{\boldsymbol{\nu}}\right) = \left[\left(\mathbf{X}^*\right)^{T} \mathbf{X}^*\right]^{-1} \left(\mathbf{X}^*\right)^{T} \hat{\boldsymbol{\nu}}^* = \left(\mathbf{X}^{T} \boldsymbol{\Omega} \mathbf{X}\right)^{-1} \mathbf{X}^T \boldsymbol{\Omega} \hat{\boldsymbol{\nu}}\]
The theoretical WV is given by
\[{\nu}_j = \left[ \frac{1}{\tau_j} \;\;\; \frac{\left(\tau_j^2 + 2\right)}{12 \tau_j} \right] \boldsymbol{\theta}.\]
Let’s simulate a process
library(wv)
library(simts)
n = 10^6
theta = c(1, 0.0000001)
model = WN(sigma2 = theta[1]) + RW(gamma2 = theta[2])
set.seed(21453)
Xt = gen_gts(n = n, model = model)
wv_Xt = wvar(Xt)
plot(wv_Xt)
Let’s add the theoretical WV:
J = floor(log2(n)) - 1
tau = 2^(1:J)
X1 = 1/tau
X2 = (tau^2 + 2)/(12*tau)
X = cbind(X1, X2)
wv_theo = X%*%theta
plot(wv_Xt, legend_position = NA)
lines(tau, wv_theo, col = "red", lwd = 2)
Let’s compute the standard GMWM estimator
Omega = diag(1/(wv_Xt$ci_high - wv_Xt$ci_low)^2)
estim_gmwm = solve(t(X)%*%Omega%*%X)%*%t(X)%*%Omega%*%wv_Xt$variance
plot(wv_Xt, legend_position = NA)
lines(tau, X%*%estim_gmwm, col = "green3", lwd = 2)