Estimating variances in time series linear regression models using empirical BLUPs and convex optimization

TL;DR

This paper introduces a novel two-stage method for estimating variance components in time series linear mixed models using empirical BLUPs and convex optimization, improving computational efficiency and theoretical understanding.

Contribution

It develops a new estimation approach for variance components in FDSLRMs, establishing theoretical existence, equivalence of estimators, and a fast, accurate algorithm for maximum likelihood estimation.

Findings

The method provides invariant, non-negative quadratic estimators applicable to any absolutely continuous distribution.

A new $ ext{O}(n)$ algorithm for (RE)MLE computation is 10^7 times more accurate and 100 times faster than existing packages.

Validated on real datasets: electricity, tourism, and cybersecurity, demonstrating practical effectiveness.

Abstract

We propose a two-stage estimation method of variance components in time series models known as FDSLRMs, whose observations can be described by a linear mixed model (LMM). We based estimating variances, fundamental quantities in a time series forecasting approach called kriging, on the empirical (plug-in) best linear unbiased predictions of unobservable random components in FDSLRM. The method, providing invariant non-negative quadratic estimators, can be used for any absolutely continuous probability distribution of time series data. As a result of applying the convex optimization and the LMM methodology, we resolved two problems $-$ theoretical existence and equivalence between least squares estimators, non-negative (M)DOOLSE, and maximum likelihood estimators, (RE)MLE, as possible starting points of our method and a practical lack of computational implementation for FDSLRM. As for computing (RE)MLE in the case of $ n $ observed time series values, we also discovered a new algorithm of order $\mathcal{O}(n)$, which at the default precision is $10^7$ times more accurate and $n^2$ times faster than the best current Python(or R)-based computational packages, namely CVXPY, CVXR, nlme, sommer and mixed. We illustrate our results on three real data sets $-$ electricity consumption, tourism and cyber security $-$ which are easily available, reproducible, sharable and modifiable in the form of interactive Jupyter notebooks.

Tables

Table 1. Table 1: First and second order moment characteristics for BLUP-NE and NE.

characteristic	NE estimators ${\breve{\sigma }}_j^2\left(X\right)$	BLUP-NE estimators $\mathring{\sigma }_j^{}{}_{}{}^2\left(X\right)$
${\Delta }_{𝛎}$	${\sigma }_0^2{\left(W^{-1}\right)}_{jj}$	$-{\sigma }_0^2{\left({𝕎}^{-1}\right)}_{jj}$
$D_{𝛎}$	$2{\left({\sigma }_j^2+{\sigma }_0^2{\left(W^{-1}\right)}_{jj}\right)}^2$	$2{\left({\sigma }_j^2-{\sigma }_0^2{\left({𝕎}^{-1}\right)}_{jj}\right)}^2$
$MS{E}_{𝛎}$	$2{\left({\sigma }_j^2+{\sigma }_0^2{\left(W^{-1}\right)}_{jj}\right)}^2+{\sigma }_0^4{\left(W^{-1}\right)}_{jj}^2$	$2{\left({\sigma }_j^2-{\sigma }_0^2{\left({𝕎}^{-1}\right)}_{jj}\right)}^2+{\sigma }_0^4{\left({𝕎}^{-1}\right)}_{jj}^2$

Table 2. Table 2: Estimation methods chosen for the second stage in EBLUP-NE.

least squares method in FDSLRM	maximum likelihood method in FDSLRM
double ordinary least squares estimators	maximum likelihood estimators
(DOOLSE) ref. stulajter_estimation_2004 , stulajter_predictions_2002	(MLE) ref. stulajter_estimation_2004 , stulajter_predictions_2002
modified (unbiased) DOOLSE	restricted (residual) MLE
(MDOOLSE) ref. stulajter_estimation_2004	(REMLE) ref. stulajter_estimation_2004
natural estimators (NE) ref. hancova_natural_2008

Table 3. Table 3: Scheme of the KKT algorithm for NN-(M)DOOLSE 𝛎 ~ + subscript ~ 𝛎 \tilde{\boldsymbol{\upnu}}_{+} in orthogonal FDSLRM.

Input: Form the matrix $G$, the vector $𝒒$ from vectors $𝒆$ and ${𝒗}_j;j=1,\dots ,l$.

For each auxiliary vector $𝒃={(b_1,b_2,\dots ,b_l)}^{\prime }\in {\{0,1\}}^l$ do: 1. Set the KKT-conditions matrix $K$: $K\leftarrow G$; For $j=1,2,\dots ,l$ do: If $b_j=0$ then $K_{0j}\leftarrow 0,K_{jj}\leftarrow -1$. 2. Calculate the auxiliary vector $𝛄$: $𝛄\leftarrow K^{-1}𝒒$. 3. Test non-negativity of $𝛄$: If $𝛄\ge 0$ then quit.

Output: Use the last $𝒃,𝛄$ to form NN-(M)DOOLSE ${\widetilde{𝛎}}_{+}$ of $𝛎$:

${\widetilde{𝛎}}_{+}\leftarrow 𝛄$;

For $j=1,2,\dots ,l$ do: If $b_j=0$ then ${\tilde{\nu }}_{+j}\leftarrow 0$.

Table 4. Table 4: Results for initial estimates and corresponding EBLUP-NE. Results for (RE)MLE are identical with NN-(M)DOOLSE.

electricity consumption - toy model 1
method	estimate $\widetilde{𝛎}$	$\Vert \widetilde{𝛎}\Vert$	EBLUP-NE $\mathring{𝛎}$	$\Vert \mathring{𝛎}\Vert$
NE	${(3.53,0.37,1.86,0.00,1.26)}^{\prime }$	$4.20$	${(3.53,0.12,1.39,0.00,0.83)}^{\prime }$	$3.89$
MLE	${(2.86,0.13,1.62,0.00,1.03)}^{\prime }$	$3.45$	${(3.53,0.05,1.42,0.00,0.84)}^{\prime }$	$3.90$
REMLE	${(3.34,0.09,1.59,0.00,0.99)}^{\prime }$	$3.83$	${(3.53,0.02,1.35,0.00,0.77)}^{\prime }$	$3.86$
electricity consumption - toy model 2
method	estimate $\widetilde{𝛎}$	$\Vert \widetilde{𝛎}\Vert$	EBLUP-NE $\mathring{𝛎}$	$\Vert \mathring{𝛎}\Vert$
NE	${(1.09,2.97,1.76,0.37,1.86)}^{\prime }$	$4.09$	${(1.09,2.79,1.59,0.24,1.69)}^{\prime }$	$3.80$
MLE	${(0.93,2.89,1.68,0.29,1.79)}^{\prime }$	$3.91$	${(1.09,2.81,1.61,0.23,1.71)}^{\prime }$	$3.83$
REMLE	${(1.09,2.87,1.67,0.28,1.77)}^{\prime }$	$3.93$	${(1.09,2.79,1.58,0.21,1.69)}^{\prime }$	$3.79$
tourism
NE	${(0.108,0.004,0.230,0.022)}^{\prime }$	$0.255$	${(0.108,0.001,0.225,0.020)}^{\prime }$	$0.250$
MLE	${(0.103,0.001,0.228,0.021)}^{\prime }$	$0.251$	${(0.108,0.000,0.225,0.020)}^{\prime }$	$0.250$
REMLE	${(0.108,0.001,0.227,0.021)}^{\prime }$	$0.253$	${(0.108,0.000,0.225,0.020)}^{\prime }$	$0.250$
cyber attacks
NE	${(0.0593,0.0255,0.0155)}^{\prime }$	$0.0664$	${(0.0593,0.0225,0.0127)}^{\prime }$	$0.0647$
MLE	${(0.0560,0.0239,0.0139)}^{\prime }$	$0.0624$	${(0.0593,0.0225,0.0125)}^{\prime }$	$0.0647$
REMLE	${(0.0593,0.0238,0.0138)}^{\prime }$	$0.0654$	${(0.0593,0.0223,0.0124)}^{\prime }$	$0.0646$

Table 5. Table 5: List of acronyms used in the paper.

acronym	explanation	specification
BLUE	best linear unbiased estimator	(3.2), p. ​3.2
BLUP	best linear unbiased predictor	(3.2), p. ​3.2
BLUP-NE	natural estimators based on BLUP	(3.1), p. ​3.1
CVXPY	Python library for convex optimization	p. ​4.2
CVXR	R version of CVXPY	p. ​4.2
DOOLSE	double ordinary least squares estimator without non-negativity constraints	(4.15), p. ​4.15
EBLUP	empirical (plug-in) BLUP	p. ​4.1
EBLUP-NE	natural estimators based on EBLUP	p. ​4.1
FDSLRM	finite discrete spectrum linear regression model	(1.1), p. ​1.1
KKT	Karush-Kuhn-Tucker	p. ​4.2
LMM	linear mixed model	p. ​2.1
LRM	linear regression model	p. ​1
MDOOLSE	modified (unbiased) DOOLSE without non-negativity constraints	(4.15), p. ​4.15
(M)DOOLSE	considering both DOOLSE and MDOOLSE	Remark 5, p. ​5
MLE	maximum likelihood estimators	(4.25), p. ​4.25
MME	Henderson’s mixed model equations	(3.2), p. ​3.2
MMEinR	R version of MATLAB function mixed	p. ​4.3
MSE	mean squared error	p. ​1
NE	natural estimators	p. ​3.1, (4.1), p. ​4.1
nlme	R package for (non) linear mixed(-effects) models	p. ​4.3
NN-DOOLSE	non-negative DOOLSE	(4.10), p. ​4.10
NN-MDOOLSE	non-negative modified DOOLSE	(4.14), p. ​4.14
NN-(M)DOOLSE	considering both NN-DOOLSE and NN-MDOOLSE	Remark 5, p. ​5
OLS	ordinary least squares	p. ​4.1
REMLE	residual (restricted) maximum likelihood estimator	(4.29), p. ​4.29
(RE)MLE	considering both MLE and REMLE	Remark 5, p. ​5
SageMath	free Python-based mathematics software	p. ​2
sommer	R package for multivariate LMMs solving MME	p. ​4.3
SciPy	Scientific Python, Python-based ecosystem of open software	p. ​2

Equations

$$X(t)=\sum _{i=1}^k{\beta }_i{f}_i(t)+\sum _{j=1}^l{Y}_j{v}_j(t)+w(t);t\in \text{mathds}T,$$

$$\begin{array}{c}\mathbf{X}=\mathrm{F}\mathbf{\beta }+\mathrm{V}\mathbf{Y}+\mathbf{w}\text{ with }E\left\{\mathbf{w}\right\}=0_n,Cov\{\mathbf{w}\}={\sigma }_0^2{\mathrm{I}}_n,\\ Cov\{\mathbf{Y},\mathbf{w}\}=0_{l\times n},Cov\{\mathbf{X}\}\equiv {\mathrm{\Sigma }}_{\text{upnu}}={\sigma }_0^2{\mathrm{I}}_n+\mathrm{VD}{\mathrm{V}}^{\prime },\end{array}$$

$$\text{identifiability}:\text{the full column rank }r(\mathrm{F}\mathrm{V})=k+l\text{ and }n>k+l.$$

$$\breve{\mathbf{Y}}={\mathrm{W}}^{-1}{\mathrm{V}}^{\prime }{\mathrm{M}}_{\mathrm{F}}\mathbf{X};\mathrm{W}={\mathrm{V}}^{\prime }{\mathrm{M}}_{\mathrm{F}}\mathrm{V}\in \text{mathds}{R}^{l\times l}.$$

$${{\mathring{\sigma }}_j}^2(\mathbf{X})=({\mathbf{Y}}_{\text{upnu}}^{\ast }{)}_j^2;{\mathbf{Y}}_{\text{upnu}}^{\ast }\text{ is the BLUP of }\mathbf{Y}\text{ based on }\mathbf{X}$$

$${{\mathring{\sigma }}_0}^2(\mathbf{X})=\frac{1}{n-k-l}[\mathbf{X}-\mathrm{F}{\mathbf{\beta }}_{\text{upnu}}^{\ast }-\mathrm{V}{\mathbf{Y}}_{\text{upnu}}^{\ast }{]}^{\prime }[\mathbf{X}-\mathrm{F}{\mathbf{\beta }}_{\text{upnu}}^{\ast }-\mathrm{V}{\mathbf{Y}}_{\text{upnu}}^{\ast }],$$

$$\left(\begin{array}{cc}{\mathrm{F}}^{\prime }\mathrm{F}& {\mathrm{F}}^{\prime }\mathrm{V}\\ {\mathrm{V}}^{\prime }\mathrm{F}& {\mathbb{G}}_{\mathrm{V}}\end{array}\right)\cdot \left(\begin{array}{c}{\mathbf{\beta }}_{\text{upnu}}^{\ast }\\ {\mathbf{Y}}_{\text{upnu}}^{\ast }\end{array}\right)=\left(\begin{array}{c}{\mathrm{F}}^{\prime }\mathbf{X}\\ {\mathrm{V}}^{\prime }\mathbf{X}\end{array}\right),$$

$${\left(\begin{array}{cc}{\mathrm{F}}^{\prime }\mathrm{F}& {\mathrm{F}}^{\prime }\mathrm{V}\\ {\mathrm{V}}^{\prime }\mathrm{F}& {\mathbb{G}}_{\mathrm{V}}\end{array}\right)}^{-1}=\left(\begin{array}{cc}({\mathrm{F}}^{\prime }\mathrm{F}{)}^{-1}& 0_{k\times l}\\ 0_{l\times k}& 0_{l\times l}\end{array}\right)+\left(\begin{array}{c}-({\mathrm{F}}^{\prime }\mathrm{F}{)}^{-1}{\mathrm{F}}^{\prime }\mathrm{V}\\ {\mathrm{I}}_l\end{array}\right){\mathbb{W}}^{-1}\left(-{\mathrm{V}}^{\prime }\mathrm{F}({\mathrm{F}}^{\prime }\mathrm{F}{)}^{-1}{\mathrm{I}}_l\right),$$

$$\left(\begin{array}{c}{\mathbf{\beta }}_{\text{upnu}}^{\ast }\\ {\mathbf{Y}}_{\text{upnu}}^{\ast }\end{array}\right)={\left(\begin{array}{cc}{\mathrm{F}}^{\prime }\mathrm{F}& {\mathrm{F}}^{\prime }\mathrm{V}\\ {\mathrm{V}}^{\prime }\mathrm{F}& {\mathbb{G}}_{\mathrm{V}}\end{array}\right)}^{-1}\cdot \left(\begin{array}{c}{\mathrm{F}}^{\prime }\mathbf{X}\\ {\mathrm{V}}^{\prime }\mathbf{X}\end{array}\right)=\left(\begin{array}{cc}\bullet & \bullet \\ -{\mathbb{W}}^{-1}{\mathrm{V}}^{\prime }\mathrm{F}({\mathrm{F}}^{\prime }\mathrm{F}{)}^{-1}& {\mathbb{W}}^{-1}\end{array}\right)\cdot \left(\begin{array}{c}{\mathrm{F}}^{\prime }\mathbf{X}\\ {\mathrm{V}}^{\prime }\mathbf{X}\end{array}\right)$$

$${\mathbf{Y}}_{\text{upnu}}^{\ast }={\mathbb{W}}^{-1}{\mathrm{V}}^{\prime }({\mathrm{I}}_n-\mathrm{F}({\mathrm{F}}^{\prime }\mathrm{F}{)}^{-1}{\mathrm{F}}^{\prime })\mathbf{X}={\mathbb{W}}^{-1}{\mathrm{V}}^{\prime }{\mathrm{M}}_{\mathrm{F}}\mathbf{X}$$

$$\mathbf{X}=\mathrm{F}\mathbf{\beta }+\mathrm{V}\mathbf{Y}+\mathbf{w},E\left\{\mathbf{w}\right\}=0_n,Cov\{\mathbf{w}\}={\sigma }_0^2{\mathrm{I}}_n,$$

$$Cov\{\mathbf{Y}\}=\mathrm{diag}\left\{{\sigma }_j^2\right\},Cov\{\mathbf{Y},\mathbf{w}\}=0_{l\times n}.$$

$${{\mathring{\sigma }}_j}^2(\mathbf{X})={({\mathbf{Y}}^{\ast })}_j^2=(\mathbb{T}\mathbf{X}{)}_j^{\prime }{(\mathbb{T}\mathbf{X})}_j={\mathbf{X}}^{\prime }{t}_j{t}_j^{\prime }\mathbf{X};j=1,\dots ,l,$$

$$\left(\begin{array}{cc}{\mathrm{F}}^{\prime }\mathrm{F}& {\mathrm{F}}^{\prime }\mathrm{VD}\\ {\mathrm{V}}^{\prime }\mathrm{F}& {\mathbb{H}}_{\mathrm{V}}\end{array}\right)\cdot \left(\begin{array}{c}{\mathbf{\beta }}_{\text{upnu}}^{\ast }\\ {\mathbf{Z}}_{\text{upnu}}^{\ast }\end{array}\right)=\left(\begin{array}{c}{\mathrm{F}}^{\prime }\mathbf{X}\\ {\mathrm{V}}^{\prime }\mathbf{X}\end{array}\right),$$

$$\mathbb{T}=\mathrm{D}{\mathbb{U}}^{-1}{\mathrm{V}}^{\prime }{\mathrm{M}}_{\mathrm{F}}\text{ with }\mathbb{U}=\mathrm{WD}+{\sigma }_0^2{\mathrm{I}}_l={\mathrm{V}}^{\prime }{\mathrm{M}}_{\mathrm{F}}\mathrm{VD}+{\sigma }_0^2{\mathrm{I}}_l.$$

$${\mathbb{H}}_{\mathrm{V}}={\mathbb{G}}_{\mathrm{V}}\mathrm{D},{\mathbb{G}}_{\mathrm{V}}^{-1}=\mathrm{D}{\mathbb{H}}_{\mathrm{V}}^{-1},\mathbb{U}=\mathbb{W}\mathrm{D},{\mathbb{W}}^{-1}=\mathrm{D}{\mathbb{U}}^{-1}.$$

$${\text{Updelta}}_{\text{upnu}}\left\{{{\mathring{\sigma }}_j}^2(\mathbf{X})\right\}=-{\sigma }_0^2({\mathbb{W}}^{-1}{)}_{jj}$$

$${\mathrm{F}}^{\prime }\mathrm{V}=0\text{ and }{\mathrm{G}}_{\mathrm{V}}={\mathrm{V}}^{\prime }\mathrm{V}=\mathrm{diag}\left\{{\Vert {\mathbf{v}}_j\Vert }^2\right\}.$$

$${\mathrm{M}}_{\mathrm{F}}\mathrm{V}=\mathrm{V},\mathbb{U}={\mathrm{G}}_{\mathrm{V}}\mathrm{D}+{\sigma }_0^2{\mathrm{I}}_l=\mathrm{diag}\left\{{\sigma }_0^2+{\sigma }_j^2{\Vert {\mathbf{v}}_j\Vert }^2\right\},$$

$${\mathbb{W}}^{-1}=\mathrm{D}{\mathbb{U}}^{-1}=\mathrm{diag}\left\{{\rho }_j/{\Vert {\mathbf{v}}_j\Vert }^2\right\},t_j={\rho }_j{\mathbf{v}}_j^{\prime }/{\Vert {\mathbf{v}}_j\Vert }^2,$$

$${{\mathring{\sigma }}_j}^2(\mathbf{X})={\rho }_j^2{\mathbf{X}}^{\prime }{\mathbf{v}}_j{\mathbf{v}}_j^{\prime }\mathbf{X}/{\Vert {\mathbf{v}}_j\Vert }^4,j=1,2,\dots ,l,$$

$${{\mathring{\sigma }}_j}^2(\mathbf{X})={\rho }_j^2{\breve{\sigma }}_j^2(\mathbf{X}),j=0,1,\dots ,l.$$

$$\breve{\text{upnu}}(\mathbf{e})=\left(\begin{array}{c}\frac{1}{n-k-l}{\mathbf{e}}^{\prime }{\mathrm{M}}_{\mathrm{V}}\mathbf{e}\\ ({\mathbf{e}}^{\prime }{\mathbf{v}}_1{)}^2/{\Vert {\mathbf{v}}_1\Vert }^4\\ ({\mathbf{e}}^{\prime }{\mathbf{v}}_2{)}^2/{\Vert {\mathbf{v}}_2\Vert }^4\\ ⋮\\ ({\mathbf{e}}^{\prime }{\mathbf{v}}_l{)}^2/{\Vert {\mathbf{v}}_l\Vert }^4\end{array}\right)=\left(\begin{array}{c}\frac{1}{n-k-l}\left({\mathbf{e}}^{\prime }\mathbf{e}-\sum _{j=1}^l({\mathbf{e}}^{\prime }{\mathbf{v}}_j{)}^2/{\Vert {\mathbf{v}}_j\Vert }^2\right)\\ ({\mathbf{e}}^{\prime }{\mathbf{v}}_1{)}^2/{\Vert {\mathbf{v}}_1\Vert }^4\\ ({\mathbf{e}}^{\prime }{\mathbf{v}}_2{)}^2/{\Vert {\mathbf{v}}_2\Vert }^4\\ ⋮\\ ({\mathbf{e}}^{\prime }{\mathbf{v}}_l{)}^2/{\Vert {\mathbf{v}}_l\Vert }^4\end{array}\right),$$

$$\displaystyle\begin{array}[]{ll}\textbf{NE}&\\ \emph{minimize}&\quad f_{0}(\boldsymbol{\upnu})=||\boldsymbol{e}\boldsymbol{e}^{\prime}-\sum\limits_{j=1}^{l}\upnu_{j}\mathrm{V}_{j}||^{2}+\left\|\mathrm{M_{V}}\boldsymbol{e}\boldsymbol{e}^{\prime}\mathrm{M_{V}}-\upnu_{0}\mathrm{V}_{0}\right\|^{2}\\[6.0pt] \emph{subject to}&\quad\boldsymbol{\upnu}=\left(\upnu_{0},\ldots,\upnu_{l}\right)^{\prime}\in[0,\infty)^{l+1}\end{array}$$

$${\mathrm{V}}_0={\mathrm{M}}_{\mathrm{F}}{\mathrm{M}}_{\mathrm{V}}\text{ and }{\mathrm{V}}_j={\mathbf{v}}_j{\mathbf{v}}_j^{\prime },j=1,\dots ,l.$$

$$\displaystyle\begin{array}[]{ll}\textbf{DOOLSE\phantom{M}}&\\[6.0pt] \emph{minimize}&f_{0}(\boldsymbol{\upnu})=\left\|\boldsymbol{e}\boldsymbol{e}^{\prime}-\mathrm{\Sigma_{\boldsymbol{\upnu}}}\right\|^{2}\\[6.0pt] \emph{subject to}&\boldsymbol{\upnu}=\left(\upnu_{0},\ldots,\upnu_{l}\right)^{\prime}\in[0,\infty)^{l+1}\end{array}$$

$$\displaystyle\begin{array}[]{ll}\textbf{MDOOLSE}&\\[6.0pt] \emph{minimize}&f_{0}(\boldsymbol{\upnu})=\left\|\boldsymbol{e}\boldsymbol{e}^{\prime}-\mathrm{M_{F}}\mathrm{\Sigma_{\boldsymbol{\upnu}}}\mathrm{M_{F}}\right\|^{2}\\[6.0pt] \emph{subject to}&\boldsymbol{\upnu}=\left(\upnu_{0},\ldots,\upnu_{l}\right)^{\prime}\in[0,\infty)^{l+1}\end{array}$$

$$\widetilde{\text{upnu}}(\mathbf{e})={\mathrm{G}}^{-1}\mathbf{q},$$

$$\displaystyle\begin{array}[]{ll}\textbf{NN-(M)DOOLSE}&\\[6.0pt] \emph{minimize}&\quad f_{0}(\boldsymbol{\upnu})=\boldsymbol{\upnu}^{\prime}\mathrm{G}\boldsymbol{\upnu}-2\boldsymbol{q}^{\prime}\boldsymbol{\upnu}\\[6.0pt] \emph{subject to}&\quad-\mathrm{I}_{l+1}\boldsymbol{\upnu}\preceq\boldsymbol{0}_{l+1}\end{array}$$

$$l_{\mathrm{m}}(\text{upnu},\mathbf{x})=\frac{1}{2}\left(\ln \det ({\mathrm{\Sigma }}_{\text{upnu}}^{-1})-{\Vert \mathbf{x}-\mathrm{F}{\mathbf{\beta }}_{\text{upnu}}^{\ast }\Vert }_{{\mathrm{\Sigma }}_{\text{upnu}}^{-1}}^2\right)$$

$$l_{\mathrm{r}}(\text{upnu},\mathbf{x})=\frac{1}{2}\left(\ln \det ({\mathrm{\Sigma }}_{\text{upnu}}^{-1})-\ln \det ({\mathrm{F}}^{\prime }{\mathrm{\Sigma }}_{\text{upnu}}^{-1}\mathrm{F})-{\Vert \mathbf{x}-\mathrm{F}{\mathbf{\beta }}_{\text{upnu}}^{\ast }\Vert }_{{\mathrm{\Sigma }}_{\text{upnu}}^{-1}}^2\right)$$

$$\displaystyle\begin{array}[]{ll}\textbf{MLE}&\\[6.0pt] \emph{maximize}&\quad f_{0}(\boldsymbol{\upnu})=\ln\det(\mathrm{\Sigma_{\boldsymbol{\upnu}}^{-1}})-\left\|\boldsymbol{e}\right\|_{\mathrm{\Sigma_{\boldsymbol{\upnu}}^{-1}}}^{2}\\[6.0pt] \emph{subject to}&\quad\boldsymbol{\upnu}=\left(\upnu_{0},\ldots,\upnu_{l}\right)^{\prime}\in[0,\infty)^{l+1}\end{array}$$

$$\displaystyle\begin{array}[]{ll}\textbf{REMLE}&\\[6.0pt] \emph{maximize}&\quad f_{0}(\boldsymbol{\upnu})=\ln\det(\mathrm{\Sigma_{\boldsymbol{\upnu}}^{-1}})-\ln\det(\mathrm{F}^{\prime}\mathrm{\Sigma_{\boldsymbol{\upnu}}^{-1}}\mathrm{F})-\left\|\boldsymbol{e}\right\|_{\mathrm{\Sigma_{\boldsymbol{\upnu}}^{-1}}}^{2}\\[6.0pt] \emph{subject to}&\quad\boldsymbol{\upnu}=\left(\upnu_{0},\ldots,\upnu_{l}\right)^{\prime}\in[0,\infty)^{l+1}\end{array}$$

Full text

1 Introduction

The need to obtain sufficiently accurate predictions for facilitating and improving decision making becomes an integral part not only of science, industry or economy but also in many other human activities. Our recent article gajdos_kriging_2017 summarizes the guiding methodology and corresponding references dealing with kriging for time series econometric forecasting as one of the advanced alternative approaches to the most popular Box-Jenkins methodology (box_time_2015 , shumway_time_2017 ).

The key idea of the time series kriging is to model the given time series data in an appropriate general class of linear regression models (LRMs) and subsequently finding the best linear unbiased predictor — the BLUP which minimizes the mean squared error (MSE) of prediction among all linear unbiased predictors, see e.g. stulajter_predictions_2002 , kedem_regression_2005 , christensen_plane_2011 , brockwell_introduction_2016 .

In the frame of kriging, we investigate theoretical features and econometric applications of a class of time series models called finite discrete spectrum linear regression models or shortly FDSLRMs. The FDSLRM class was introduced in 2002-2003 by Štulajter (stulajter_predictions_2002 , stulajter_mse_2003 ) as a direct extension of classical (ordinary) regression models (see e.g. christensen_plane_2011 , hyndman_forecasting:_2018 ).

The FDSLRM has mean values (trend) given by linear regression and random components (error terms) are represented as a sum of a linear combination of uncorrelated zero-mean random variables and white noise which can be interpreted in terms of the finite discrete spectrum priestley_spectral_2004 .

Formally the FDSLRM can be presented as

[TABLE]

where

$\mathds{T}$ representing the time domain is a countable subset of the real line $\mathds{R}$,

$k$ and $l$ are some fixed non-negative integers, i.e. $k,l\in\mathds{N}_{0}$,

$\boldsymbol{\beta}=(\beta_{1},\beta_{2},\ldots,\beta_{k})^{\prime}\in\mathds{R}^{k}$ is a vector of regression parameters,

$\boldsymbol{Y}=(Y_{1},Y_{2},\ldots,Y_{l})^{\prime}$ is an unobservable $l\times 1$ random vector with zero mean vector $E\left\{\boldsymbol{Y}\right\}=\boldsymbol{0}_{l}$, and with an $l\times l$ diagonal covariance matrix ${Cov}\{\boldsymbol{Y}\}=\mathrm{diag}\left\{\sigma_{j}^{2}\right\}$, where $\sigma_{j}^{2}\in\mathds{R}_{+}$ are non-negative real numbers,

$f_{i}(.);\,i=1,2,\ldots,k$ and $v_{j}(.);\,j=1,2,\ldots,l$ are real functions defined on $\mathds{R}$,

$w(.)$ stands for white noise uncorrelated with $\boldsymbol{Y}$ and having a positive dispersion $D\left\{w(t)\right\}=\sigma_{0}^{2}\in\mathds{R}_{++}$.

Typically, due to the nature of time series data collection, the most frequently considered time domain $\mathds{T}$ is the set of natural numbers $\mathds{N}=\{1,2,\ldots\}$. FDSLRM variance parameters, which are fundamental quantities in kriging, are commonly described by one vector $\boldsymbol{\upnu}\equiv(\upnu_{0},\upnu_{1},\ldots,\upnu_{l})^{\prime}=(\sigma_{0}^{2},\sigma_{1}^{2},\ldots,\sigma_{l}^{2})^{\prime}$, an element of the parametric space $\Upupsilon=(0,\infty)\times[0,\infty)^{l}$ or $\Upupsilon=\mathds{R}_{++}\times\mathds{R}^{l}_{+}$ for short.

The principal goal of our paper is to present an alternative estimation method for FDSLRM variances, one of the time series kriging steps gajdos_kriging_2017 . Our approach represents an improved two-stage modification of the previously developed method of natural estimators (NE) hancova_natural_2008 which is now based on the idea of empirical (plug-in) BLUPs. We will refer to our new method as EBLUP-NE for short.

The structure of the paper is organized as follows. Section 2 includes several important notes about theoretical methods and computational tools for FDSLRM kriging offered by closely connected mathematical branches and computational research. This background became the strong basis for our research.

Section 3 describes the first stage of our EBLUP-NE method. It develops the definition and computational form of EBLUP-NE at known variance parameters $\boldsymbol{\upnu}$ using two special matrices known in linear algebra as the Schur complement and Gram matrix zhang_schur_2005 . The method requires very minimal distributional assumptions on a finite FDLSRM observation $\boldsymbol{X}=(X(1),\ldots,X(n))^{\prime},$ $n\in\mathds{N}$. It must have only an absolutely continuous probability distribution with respect to some $\sigma$-finite measure. In section 3 we also derive the basic statistical properties of EBLUP-NE at known $\boldsymbol{\upnu}$.

Section 4 is devoted to the second stage of EBLUP-NE method which consists in replacing true, in practice unknown variance parameters $\boldsymbol{\upnu}$ by another, appropriate FDSLRM estimates of $\boldsymbol{\upnu}$ obtained by maximum likelihood or least squares gajdos_kriging_2017 . Since the theory and computational aspects of these estimation procedures in FDSLRM fitting based on the projection theory in Hilbert spaces are more than 10 years old (stulajter_predictions_2002 , stulajter_estimation_2004 , hancova_natural_2008 ), we revisit and update them in the light of recent advances in closely related linear mixed modeling and convex optimization.

In the following fifth section, we illustrate theoretical results of the paper and the performance of EBLUP-NE on three real time series data sets — electricity consumption, tourism and cyber security. Section 6 presents the conclusions of the paper. For the sake of paper readability, we moved very technical details and proofs to the Appendix. In the Appendix we also report a list of acronyms and abbreviations used in the paper (tab. 5).

Since the paper is aimed at statisticians, analysts, econometricians and data scientists applying time series and forecasting methods, our notation is standard for time series analysis and prediction using linear regression models (stulajter_predictions_2002 , kedem_regression_2005 , brockwell_time_2009 ). Sets like $\mathds{R},\mathds{R}_{+},\mathds{R}_{++}$ are labeled as it is common in convex optimization boyd_convex_2009 .

2 Theoretical and computational bases for FDSLRM kriging

Let us recall some important notes about time series FDSLRM kriging and its connection to other math branches and current computational technology.

From a modeling point of view, if we consider any time series model $X(.)$ in the additive form $X(t)=m(t)+\varepsilon(t);t\in\mathds{T}$ whose mean value function $m(.)$ is a real function on $\mathds{R}$ expressible by a functional series (e.g. Taylor or Fourier series) and error term $\varepsilon(.)$ is a mean-zero stationary process, then according to the spectral representation theory of time series (priestley_spectral_2004 ,percival_spectral_2009 ) $X(.)$ can be approximated arbitrarily closely by FDSLRM. Therefore, from a practice perspective, FDSLRMs can be potentially applied in many practical situations.

In econometric applications of FDSLRM kriging, we almost always have only one realization of time series $X(.)$ and we do not know the mean-value parameters $\boldsymbol{\beta}\in\mathds{R}^{k}$ nor the variance parameters $\boldsymbol{\upnu}\in\Upupsilon$. From the view of time series analysis, the problem of estimating variances $\boldsymbol{\upnu}$ belongs to the class of problems concerning covariance matrix estimation with one realization wu_covariance_2012 .

In addition, any finite FDSLRM observation $\boldsymbol{X}$ satisfies a special type of linear mixed model (LMM) of the form

[TABLE]

where design matrices $\mathrm{F}=\{F_{ti}\}=\{f_{i}(t)\},\mathrm{V}=\{V_{tj}\}=\{v_{j}(t)\}\text{ for }t=1,2,\ldots,n$; $i=1,2,\ldots,k$; $j=1,2\ldots,l$ and random vector $\boldsymbol{w}=(w(1),w(2),\ldots,w(n))^{\prime}$ is a finite $n$-dimensional white noise observation.

The fundamental property (2.1) of FDSLRM allows us to apply convenient LMM mathematical techniques for FDSLRM fitting and forecasting (pinheiro_mixed-effects_2009 , christensen_plane_2011 , witkovsky_estimation_2012 , demidenko_mixed_2013 , rao_small_2015 , covarrubias-pazaran_genome-assisted_2016 , singer_graphical_2017 ) together with up-to-date LMM software packages — nlme, sommer written in R (pinheiro_nlme:_2018 , covarrubias-pazaran_sommer:_2019 ), MATLAB function mixed witkovsky_mixed_2018 or SAS package Proc MIXED sasinstituteinc._sas/stat_2018 . As the first practical consequence of the LMM theory, to have identifiable FDSRLMs under assumptions of multivariate normal distribution (demidenko_mixed_2013, , sec. 3.2, thm 11), in the paper, we will assume for the block matrix $(\mathrm{F}\;\mathrm{V})$ and for the number of data $n$ in $\boldsymbol{x}$, a realization of $\boldsymbol{X}$ in real situations, the following sufficient condition

[TABLE]

On the other hand, the standard maximum likelihood or least squares estimates of $\boldsymbol{\upnu}$ in FDSLRM are optimization problems. Our motivation to explore estimating FDSLRM variances in the frame of convex optimization lays on the fact that its mathematical tools and efficient, very reliable computational interior-point methods became important or fundamental tools in many other branches of mathematics boyd_convex_2009 like the design of experiments, high-dimensional data, machine learning or data mining.

Inspired by essential and recent works on convex optimization (boyd_convex_2009 , bertsekas_convex_2009 , koenker_convex_2014 , cornuejols_optimization_2018 , agrawal_rewriting_2018 ) we prove new theoretical relations among existing FDSLRM estimators in the extended parametric space $\Upupsilon=\mathds{R}^{l+1}_{+}$ for the orthogonal version of FDSLRM, most used in real time series applications. Moreover, we show how to apply the latest convex optimization packages CVXPY and CVXR based on disciplined convex programming, written in Python diamond_cvxpy:_2016 and R fu_cvxr:_2019 , for estimating $\boldsymbol{\upnu}$. Finally, we also focus on the development of a new fast and accurate computational optimization algorithm for estimating variances in FDSLRM.

As for computational technology, our time series calculations are carried out using free open source software based on the R statistical language and packages (rdevelopmentcoreteam_r:_2019 , mcleod_time_2012 ), the Scientific Python with the SciPy ecosystem (jones_scipy:_2001 , oliphant_python_2007 ) and free Python-based mathematics software SageMath (stein_sage_2019 , hogben_sage_2013 , zimmermann_computational_2018 ), an open source alternative to the well-known commercial computer algebra systems Mathematica or Maple. The simultaneous use of R, SciPy and SageMath provides us valuable cross-checking of our computational results.

All our computational algorithms and results are easily readable, sharable, reproducible, and modifiable thanks to open-source Jupyter technology kluyver_jupyter_2016 . In particular, we present our results in the form of Jupyter notebooks, dynamic HTLM documents integrating prose, code and results similarly as Mathematica notebooks, stored in our free available GitHub repository (gajdos_fdslrm_2019 , adar_coding_2014 ).

The Jupyter notebooks, detailed records of our computing with explaining narratives, can be seen or studied as static HTML pages via Jupyter nbviewer (https://nbviewer.jupyter.org/, kluyver_jupyter_2016 ) or interactively as live HTML documents using Binder (https://mybinder.org/, projectjupyter_binder_2018 ) where the code is executable. The way and presentation of our computing with real data are inspired by works brieulle_computing_2019 , weiss_scientific_2017 .

3 The BLUP estimation method for $\boldsymbol{\upnu}$

3.1 Definition and computational form of estimators

In the paper hancova_natural_2008 , we proposed the method of always non-negative natural estimators (NE) of FDSLRM variances $\boldsymbol{\upnu}$. The main idea behind NE came from the fact that $\sigma_{j}^{2}={Cov}\{Y_{j}\}=E\left\{Y_{j}^{2}\right\};j=1,\ldots,l$. Therefore, if random vector $\boldsymbol{Y}$ in (1.1) was known, the natural estimate of $\sigma_{j}^{2}$ would be just $Y_{j}^{2}$. This initial consideration was identical with Rao’s estimates known as MINQUE ((rao_estimation_1988, , sec. 5.1), (christensen_plane_2011, , sec. 12.7)).

To get sufficiently simple, explicit analytic expressions available for further theoretical study, we predicted the random vector $\boldsymbol{Y}$ by the ordinary least squares method leading to the following linear predictor of $\boldsymbol{Y}$ based on $\boldsymbol{X}$

[TABLE]

Matrix $\mathrm{M_{F}}=(\mathrm{I}_{n}-\mathrm{F}(\mathrm{F}^{\prime}\mathrm{F})^{-1}\mathrm{F}^{\prime})\in\mathds{R}^{n\times n}$ represents (hancova_natural_2008 ) the orthogonal projector onto the orthogonal complement of the column space of $\mathrm{F}$ and matrix $\mathrm{W}\succ 0$ ($\succ 0$ means positive definiteness) is known in linear algebra and statistics (zhang_schur_2005, , chap. 6) as the Schur complement of $\mathrm{F}^{\prime}\mathrm{F}$ in the block matrix $\left(\begin{smallmatrix}\mathrm{F}^{\prime}\mathrm{F}&\mathrm{F}^{\prime}\mathrm{V}\\ \mathrm{V}^{\prime}\mathrm{F}&\mathrm{V}^{\prime}\mathrm{V}\\ \end{smallmatrix}\right)$.

However, from perspective of the general theory of the best linear unbiased prediction in LMM (see e.g. (christensen_plane_2011, , chap. 12)), we could consider natural estimators of $\boldsymbol{\upnu}$ in FDSLRM based on the BLUP $\boldsymbol{Y}_{\!\!\boldsymbol{\upnu}}^{*}$ of $\boldsymbol{Y}$. Then intuitively, it seems reasonable to assume that the BLUP as the unbiased linear predictor with minimal MSE could lead to better estimators of $\boldsymbol{\upnu}$. This intuition will be confirmed by our following derivations.

The previous heuristic consideration motivates the following definition. Let us consider FDSLRM (1.1) and its observation (2.1). Then estimators in the form

[TABLE]

are called natural estimators of $\boldsymbol{\upnu}$ based on the BLUP of $\boldsymbol{Y}$ or BLUP-NE for short.

Remark 1

Following the analogy with original NE in hancova_natural_2008 , for the first component $\upnu_{0}=\sigma_{0}^{2}$ of $\boldsymbol{\upnu}$ we could take as a BLUP based natural estimator the sum of squares of white noise residuals $\boldsymbol{w}_{\boldsymbol{\upnu}}^{*}=\boldsymbol{X}-\mathrm{F}\boldsymbol{\beta}_{\boldsymbol{\upnu}}^{*}-\mathrm{V}\boldsymbol{Y}_{\!\!\boldsymbol{\upnu}}^{*}$ divided by the number of degrees of freedom in FDSLRM

[TABLE]

where $\boldsymbol{\beta}_{\boldsymbol{\upnu}}^{*}$ is the best linear unbiased estimator (BLUE) of $\boldsymbol{\beta}$. In the LMM framework, FDSLRM white noise residuals $\boldsymbol{w}_{\boldsymbol{\upnu}}^{*}$ are also called conditional residuals singer_graphical_2017 .

But now our intuition fails. The computational form and properties of $\boldsymbol{Y}_{\!\!\boldsymbol{\upnu}}^{*}$ depend on all variance components $\boldsymbol{\upnu}$. That dependency also appears in the proposed estimator ${\mathring{{\sigma}}_{0}}^{2}(\boldsymbol{X})$ which causes its biasedness even in a simpler special case of FDSLRM (orthogonality condition, see (3.9)). Therefore, further we do not pay any special attention to its form and properties. As a natural estimator of $\sigma_{0}^{2}$ we will keep the original unbiased invariant quadratic NE $\breve{\sigma}_{0}^{2}(\boldsymbol{X})$ of $\sigma_{0}^{2}$, eq. (2.2) in hancova_natural_2008 .

The BLUP $\boldsymbol{Y}_{\!\!\boldsymbol{\upnu}}^{*}$ of $\boldsymbol{Y}$ together with the BLUE $\boldsymbol{\beta}_{\boldsymbol{\upnu}}^{*}$ of $\boldsymbol{\beta}$ can be obtained from celebrated Henderson’s mixed model equations, MME for short (for MME in the current framework of LMM see e.g. (christensen_plane_2011, , sec. 12.3) or (witkovsky_estimation_2012, , sec. 2)). MME have the following form in the case of FDSLRM

[TABLE]

where $\mathbb{G}_{\mathrm{V}}=\mathrm{G}_{\mathrm{V}}+\sigma_{0}^{2}\mathrm{D}^{-1}=\mathrm{V}^{\prime}\mathrm{V}+\sigma_{0}^{2}\mathrm{D}^{-1}$ and $\mathrm{G}_{\mathrm{V}}=\mathrm{V}^{\prime}\mathrm{V}\in\mathds{R}^{l\times l}$ is the Gram matrix ((boyd_introduction_2018, , sec. 10.1)) for columns of design matrix $\mathrm{V}$.

Under the assumption of the full column rank $r=k+l$ of matrix $(\mathrm{F}\;\mathrm{V})$, when FDSRLM is identifiable, and using the so-called Banachiewicz inversion formula for the inverse of a $2\times 2$ partitioned (block) matrix (see e.g. (zhang_schur_2005, , sec. 6.0.2) or (hancova_natural_2008, , sec. 2.1)), we can easily prove the existence and the following form of the inverse to the block matrix in (3.2)

[TABLE]

where $\mathbb{W}=\mathrm{W}+\sigma_{0}^{2}\mathrm{D}^{-1}=\mathrm{V}^{\prime}\mathrm{M_{F}}\mathrm{V}+\sigma_{0}^{2}\mathrm{D}^{-1}\succ 0$ is again the Schur complement of $\mathrm{F}^{\prime}\mathrm{F}$ but now in the block matrix $\left(\begin{smallmatrix}\mathrm{F}^{\prime}\mathrm{F}&\mathrm{F}^{\prime}\mathrm{V}\\ \mathrm{V}^{\prime}\mathrm{F}&\mathbb{G}_{\mathrm{V}}\\ \end{smallmatrix}\right)$. We also call matrix $\mathbb{W}$ more specifically as the extended Schur complement determined by variances $\boldsymbol{\upnu}$.

Substituting the last result for the inverse into MME (3.2) and rearranging, we get for $\boldsymbol{Y}_{\!\!\boldsymbol{\upnu}}^{*}$ (symbols $\bullet$ denote blocks not needed in deriving)

[TABLE]

[TABLE]

Denoting matrix $\mathbb{W}^{-1}\mathrm{V}^{\prime}\mathrm{M_{F}}\in\mathds{R}^{l\times n}$ as $\mathbb{T}$, we just derived the following computational form of BLUP-NE estimators.

Proposition 1 (the computational form of BLUP-NE)

Let us consider the following LMM model for FDSLRM observation $\boldsymbol{X}$

[TABLE]

Then BLUP-NE estimators $\mathring{{\sigma}}_{1}^{2},\ldots,\mathring{{\sigma}}_{l}^{2}$ of parameters $\sigma_{1}^{2},\ldots,\sigma_{l}^{2}$ are given by

[TABLE]

*where $\mathbb{t}_{j}=(\mathbb{T}_{j1},\mathbb{T}_{j2},\ldots,\mathbb{T}_{jn})^{\prime}$ are rows of matrix $\mathbb{T}$ and $\mathbb{T}=\mathbb{W}^{-1}\mathrm{V}^{\prime}\mathrm{M_{F}}$ with $\mathbb{W}=\mathrm{diag}\left\{\sigma_{0}^{2}/\sigma_{j}^{2}\right\}+\mathrm{V}^{\prime}\mathrm{M_{F}}\mathrm{V}$.

Remark 2

It is straightforward to extend our computational form to the case when some variances $\sigma_{j}^{2}$ are zero or in other words when matrix $\mathrm{D}$ is singular. For a singular $\mathrm{D}$, MME (3.2) have the alternative form witkovsky_estimation_2012

[TABLE]

where $\mathbb{H}_{\mathrm{V}}=\mathrm{G}_{\mathrm{V}}\mathrm{D}+\sigma_{0}^{2}\mathrm{I}_{l}=\mathrm{V}^{\prime}\mathrm{V}\mathrm{D}+\sigma_{0}^{2}\mathrm{I}_{l}$ and $\boldsymbol{Y}^{*}=\mathrm{D}\boldsymbol{Z}^{*}$.

Applying the same argument with the Banachiewicz inversion formula, we get the second version of the computational form of $\mathbb{T}$ determining $\boldsymbol{Y}_{\!\!\boldsymbol{\upnu}}^{*}$

[TABLE]

As opposed to $\mathbb{W},\mathbb{G}_{\mathrm{V}}$, matrices $\mathbb{U},\mathbb{H}_{\mathrm{V}}$ always exist under identifiability assumptions of our FDSLRM. Simultaneously both of $\mathbb{U},\mathbb{H}_{\mathrm{V}}$ are also always invertible, since both of $\det\mathbb{U},\det\mathbb{H}_{\mathrm{V}}$ are nonzero for any $\mathrm{D}$ in our FDSLRM as it is shown in the Appendix. In the case of a nonsingular $\mathrm{D}$ the mentioned matrices are connected via following relationships

[TABLE]

The version (3.5) of MME is also preferred in numerical calculations witkovsky_estimation_2012 , since it can handle not only a singular $\mathrm{D}$ but also a very ill-conditioned $\mathrm{D}$ appearing when $\boldsymbol{\upnu}$ has very small positive components $\sigma_{j}^{2}$.

It is worth to mention that originally MME were derived under the normality assumptions (see the original work henderson_estimation_1959 or witkovsky_estimation_2012 ), but from the viewpoint of least squares both versions (3.2), (3.5) of MME describe the BLUP $\boldsymbol{Y}_{\!\!\boldsymbol{\upnu}}^{*}$ and BLUE $\boldsymbol{\beta}_{\boldsymbol{\upnu}}^{*}$ ((christensen_plane_2011, , sec. 12.3)) with no need to restrict distributions of $\boldsymbol{Y}$ and $\boldsymbol{w}$ to be normal.

3.2 Statistical properties at known variance parameters

The derivation of theoretical properties of BLUP-NE estimates under the assumption of known $\boldsymbol{\upnu}$, regarding the first and second order moment characteristics, can be significantly facilitated by the following lemma describing the properties of matrix $\mathbb{T}$ determining the BLUP $\boldsymbol{Y}_{\!\!\boldsymbol{\upnu}}^{*}$ in (3.3).

Its proof can be accomplished in a similar way as the proof of Lemma 3.1 in hancova_natural_2008 by a direct routine computation employing formulas (3.7), $\mathbb{T}=\mathbb{W}^{-1}\mathrm{V}^{\prime}\mathrm{M_{F}}$, $\mathbb{W}=\mathrm{V}^{\prime}\mathrm{M_{F}}\mathrm{V}+\sigma_{0}^{2}\mathrm{D}^{-1}$, properties of orthogonal projectors (orthogonality, idempotency, symmetry) and Schur complements (symmetry, positive definiteness).

Lemma 1 (basic properties of $\mathbb{T}$)

1. $(0)$

$\mathbb{T}\mathbb{T}^{\prime}=\mathbb{W}^{-1}(\mathrm{I}_{l}-\sigma_{0}^{2}\mathrm{D}^{-1}\mathbb{W}^{-1})=\mathrm{D}\mathbb{U}^{-1}(\mathrm{I}_{l}-\sigma_{0}^{2}\mathbb{U}^{-1})$** 2. $(0)$

$\mathbb{T}\mathrm{F}=\boldsymbol{0}_{l\times n}$* and $\mathbb{T}\mathrm{V}=\mathrm{I}_{l}-\sigma_{0}^{2}\mathrm{D}^{-1}\mathbb{W}^{-1}=\mathrm{I}_{l}-\sigma_{0}^{2}\mathbb{U}^{-1}$* 3. $(0)$

$\mathbb{T}\mathrm{\Sigma_{\boldsymbol{\upnu}}}\mathbb{T}^{\prime}=\mathrm{D}-\sigma_{0}^{2}\mathbb{W}^{-1}=\mathrm{D}(\mathrm{I}_{l}-\sigma_{0}^{2}\mathbb{U}^{-1})$**

The computational forms (3.4) of BLUP-NE are also quadratic forms of $\boldsymbol{X}$. In addition, result (2) $\mathbb{T}\mathrm{F}=\boldsymbol{0}_{l\times n}$ implies that $\mathbb{t}_{j}\mathrm{F}=\boldsymbol{0}$. Such a condition leads to the conclusion that BLUP natural estimators ${\mathring{{\sigma}}_{j}}^{2}(\boldsymbol{X})$ are translation invariant or shortly invariant quadratic estimators (stulajter_predictions_2002, , sec. 1.5) or (hancova_natural_2008, , sec. 3.1). The following theorem summarizes theoretical properties of BLUP-NE.

Theorem 3.1 (statistical properties of ${\mathring{{\sigma}}_{j}}^{2}(\boldsymbol{X})$)

Natural estimators ${\mathring{{\sigma}}_{j}}^{2}(\boldsymbol{X});j=1,2,\ldots,l$ of $\boldsymbol{\upnu}$ based on the BLUP $\boldsymbol{Y}_{\!\!\boldsymbol{\upnu}}^{*}$ of $\boldsymbol{Y}$ are invariant quadratic estimators having the following properties

1. $(0)$

$E_{\boldsymbol{\upnu}}\left\{{\mathring{{\sigma}}_{j}}^{2}(\boldsymbol{X})\right\}=\sigma_{j}^{2}-\sigma_{0}^{2}(\mathbb{W}^{-1})_{jj};$* $j=1,2,\ldots,l,$*

If $\boldsymbol{X}\sim\mathcal{N}_{n}(\mathrm{F}\boldsymbol{\beta},\mathrm{\Sigma_{\boldsymbol{\upnu}}})$, then 2. $(0)$

$D_{\boldsymbol{\upnu}}\{{\mathring{{\sigma}}_{j}}^{2}(\boldsymbol{X})\}=2(\sigma_{j}^{2}-\sigma_{0}^{2}(\mathbb{W}^{-1})_{jj})^{2};$* $j=1,2,\ldots,l,$* 3. $(0)$

${Cov}_{\boldsymbol{\upnu}}\left\{{\mathring{{\sigma}}_{i}}^{2}(\boldsymbol{X}),{\mathring{{\sigma}}_{j}}^{2}(\boldsymbol{X})\right\}=2(\sigma_{0}^{2}(\mathbb{W}^{-1})_{ij})^{2};$* $i,j=1,2,\ldots,l,$ $i\neq j$,* 4. $(0)$

${MSE}_{\boldsymbol{\upnu}}\{{\mathring{{\sigma}}_{j}}^{2}(\boldsymbol{X})\}=2{(\sigma_{j}^{2}-\sigma_{0}^{2}(\mathbb{W}^{-1})_{jj})}^{2}+\sigma_{0}^{4}(\mathbb{W}^{-1})_{jj}^{2};$* $j=1,2,\ldots,l.$*

Proof

See the Appendix.

Remark 3

If some components in $\boldsymbol{\upnu}$ are zero, there is a need to use expression $\mathrm{D}\mathbb{U}^{-1}$ instead of $\mathbb{W}^{-1}$. The first property, biasedness of ${\mathring{{\sigma}}_{j}}^{2}(\boldsymbol{X}),j=1,2,\ldots,l$, is in full accordance with the Ghosh theorem (ghosh_nonexistence_1996 , gajdos_kriging_2017 ) about the incompatibility between simultaneous non-negativity and unbiasedness of estimators for variances of random components $\boldsymbol{Y}$ in LMM. As for bias, defined as $\Updelta_{\boldsymbol{\upnu}}\left\{{\mathring{{\sigma}}_{j}}^{2}(\boldsymbol{X})\right\}=E_{\boldsymbol{\upnu}}\left\{{\mathring{{\sigma}}_{j}}^{2}(\boldsymbol{X})\right\}-\sigma_{j}^{2}$, we have the following expression given by the extended Schur complement $\mathbb{W}$

[TABLE]

Now we can theoretically compare quality of the BLUP natural estimators with original natural estimators (NE) proposed in our previous paper hancova_natural_2008 . As we can see from summarizing tab. 1, their quality is determined by Schur complements $\mathrm{W}$ and $\mathbb{W}$ where both of them are symmetric and positive definite matrices.

Using elementary properties of the Löwner partial ordering (puntanen_formulas_2013, , chap. 24), we have immediately relations $\mathbb{W}>\mathrm{W},\mathrm{W}^{-1}>\mathbb{W}^{-1},(\mathrm{W})_{jj}>(\mathbb{W}^{-1})_{jj}$ which imply a smaller absolute value of bias, dispersion and MSE for BLUP-NE ${\mathring{{\sigma}}_{j}}^{2}(\boldsymbol{X})$ in comparison with original NE $\breve{\sigma}_{j}^{2}(\boldsymbol{X})$. This conclusion directly supports our heuristic idea leading to the BLUP-NE definition (3.1).

Orthogonal FDSLRM

In real time series analysis using FDSLRMs, the fundamental modeling procedure is based on spectral analysis of time series (priestley_spectral_2004 , brockwell_time_2009 ). We identify the significant Fourier frequencies by periodogram (stulajter_estimation_2004 , gajdos_kriging_2017 ) which restrict the form of design matrices $\mathrm{F}$ an $\mathrm{V}$ leading to the orthogonality condition for FDSLRM stulajter_estimation_2004

[TABLE]

where $\boldsymbol{v}_{j},j=1,2\ldots,l$ is j-th column of $\mathrm{V}$. Such a FDSLRM, satisfying condition (3.9), is called orthogonal stulajter_estimation_2004 . So in practice, we mainly work with the orthogonal version of FDSLRM.

Under the condition of orthogonality (3.9), we can write for matrices $\mathrm{M_{F}}$, $\mathbb{U}$, $\mathbb{W}$, $\mathbb{T}$ and BLUP-NE ${\mathring{{\sigma}}_{j}}^{2}(\boldsymbol{X})$

[TABLE]

where we introduced $\rho_{j}\equiv\frac{\sigma_{j}^{2}\left\|\boldsymbol{v}_{j}\right\|^{2}}{\sigma_{0}^{2}+\sigma_{j}^{2}\left\|\boldsymbol{v}_{j}\right\|^{2}}\in\mathds{R},\;0\leqq\rho_{j}<1,j=1,2,\ldots,l$.

Applying these results, we get the following direct corollary of Theorem 3.1 for any orthogonal FDSLRM.

Corollary 1 (properties of ${\mathring{{\sigma}}_{j}}^{2}(\boldsymbol{X})$ in orthogonal FDSLRM)

In an orthogonal FDSLRM natural estimators ${\mathring{{\sigma}}_{j}}^{2}$ based on the BLUP $\boldsymbol{Y}_{\!\!\boldsymbol{\upnu}}^{*}$ have the following properties

1. $(0)$

$E_{\boldsymbol{\upnu}}\left\{{\mathring{{\sigma}}_{j}}^{2}(\boldsymbol{X})\right\}=\rho_{j}\sigma_{j}^{2}$; $j=1,2,\ldots,l,$ where $\rho_{j}=\dfrac{\sigma_{j}^{2}\left\|\boldsymbol{v}_{j}\right\|^{2}}{\sigma_{0}^{2}+\sigma_{j}^{2}\left\|\boldsymbol{v}_{j}\right\|^{2}}\in\mathds{R}_{+}$.

If $\boldsymbol{X}\sim\mathcal{N}_{n}(\mathrm{F}\boldsymbol{\beta},\mathrm{\Sigma_{\boldsymbol{\upnu}}})$, then 2. $(0)$

$D_{\boldsymbol{\upnu}}\{{\mathring{{\sigma}}_{j}}^{2}(\boldsymbol{X})\}=2\rho_{j}^{2}\sigma^{4}_{j}$; $j=1,2,\ldots,l,$ 3. $(0)$

${Cov}_{\boldsymbol{\upnu}}\left\{{\mathring{{\sigma}}_{i}}^{2}(\boldsymbol{X}),{\mathring{{\sigma}}_{j}}^{2}(\boldsymbol{X})\right\}=0$; $i,j=1,2,\ldots,l,$ $i\neq j$, 4. $(0)$

${MSE}_{\boldsymbol{\upnu}}\{{\mathring{{\sigma}}_{j}}^{2}(\boldsymbol{X})\}=[2\rho_{j}^{2}+(1-\rho_{j})^{2}]\sigma^{4}_{j}$; $j=1,2,\ldots,l.$

Finally, if we look at original natural estimators $\breve{\sigma}_{j}^{2}(\boldsymbol{X})$ from our paper hancova_natural_2008 , then in orthogonal FDSLRM it can be easily shown that ${\mathring{{\sigma}}_{j}}^{2}(\boldsymbol{X})=\rho_{j}^{2}\breve{\sigma}_{j}^{2}(\boldsymbol{X}),j=1,2,\ldots,l$. If we introduce $\rho_{0}=1$, then we obtain the complete orthogonal version of (3.4) for computing BLUP-NE ${\mathring{{\sigma}}_{j}}^{2}(\boldsymbol{X})$:

[TABLE]

4 Mathematical and computational tools of convex optimization

4.1 Empirical BLUPs in the BLUP-NE estimation method

The computational form of our BLUP natural estimators, their first and second order moment characteristics were derived at known variances $\boldsymbol{\upnu}$. But if we really knew variances $\boldsymbol{\upnu}$, we would not need to estimate them. The practical estimation procedure for $\boldsymbol{\upnu}$ cannot depend on the parameters $\boldsymbol{\upnu}$ which are to be estimated.

However, such situation is not rare at all in LMMs, also in the FDSLRM theory, when we consider e.g. double weighted least squares estimators (DOWELSE) or maximum likelihood estimators (MLE) for $\boldsymbol{\upnu}$. In FDSRLM fitting these estimators are computed by iterative numerical procedures using projections in Hilbert spaces ((stulajter_predictions_2002, , sec. 3.4) or stulajter_estimation_2004 ). Generally, in the LMM framework, similar iterative numerical estimation procedures are described in (searle_variance_2009, , chap. 8), witkovsky_estimation_2012 .

In our case, the simplest reasonable solution of the mentioned situation is to use an empirical version of BLUP (EBLUP) based on the computational form of BLUP-NE (3.4), or (3.10) in an orthogonal case, as it is usual in the general theory of empirical BLUPs in LMMs ((stulajter_predictions_2002, , chap. 5), witkovsky_estimation_2012 , (rao_small_2015, , chap. 5)). In particular, this step means replacing unknown true parameters $\boldsymbol{\upnu}$ with other ”initial” values or estimates of $\boldsymbol{\upnu}$ in FDSLRM. In the light of iterative approaches, our EBLUP-NE can be viewed as a two-stage iterative method with one step in the iteration.

If we think again heuristically as during the formulation of the BLUP-NE definition (3.1), it seems reasonable to expect as good starting values of $\boldsymbol{\upnu}$ any, in some sense optimal estimates of $\boldsymbol{\upnu}$, e.g. estimates obtained by maximum likelihood or least squares in FDSLRM.

On the top of that, from the theoretical perspective we should prefer estimates with the simplest explicit form and under less restrictive assumptions on the structure or distribution of FDSLRM. On the other hand, from the practical point of view, it will be sufficient to have at least initial estimates which can be obtained by reliable and time efficient computational methods.

Therefore, in the following sections we investigate, theoretically and practically with real data sets, five estimation methods providing initial estimates for $\boldsymbol{\upnu}$ under different assumptions on the structure and distribution of FDSLRM observation $\boldsymbol{X}$. Our candidates for the second stage of EBLUP-NE are summarized in tab. 2.

Since theoretical and computational aspects of these estimation procedures in FDSLRM fitting based on the projection theory in Hilbert spaces are more than 10 years old (stulajter_predictions_2002 , stulajter_estimation_2004 , hancova_natural_2008 ), we revisit and update them in the light of the current theoretical and computational tools of convex optimization.

According to agrawal_rewriting_2018 , boyd_convex_2009 , to formulate any mathematical problem as an optimization problem, we need to identify three attributes of the problem: optimization variable, constraints that the variable must satisfy and the objective function depending on the variable whose optimal value we want to achieve.

In all our estimation methods, optimization variable is $\boldsymbol{\upnu}\in\mathds{R}^{l+1}$ satisfying constraints $\Upupsilon=\mathds{R}_{++}\times\mathds{R}^{l}_{+}$. In order to avoid any incompatibility problems in the convex optimization theoretical or computational framework, we extend constraints $\Upupsilon$ into the form of standard nonnegativity constraints $\Upupsilon^{*}=[0,\infty]^{l+1}$ or $\boldsymbol{\upnu}\succeq 0$ using generalized inequality.

Finally, we also consider using present theoretical and computational tools used in linear mixed modeling, since in the LMM framework, the problem of estimating variances has a long and rich history with many essential reference works (e.g. rao_estimation_1988 , searle_variance_2009 , christensen_plane_2011 , demidenko_mixed_2013 , rao_small_2015 , see also a review paper witkovsky_estimation_2012 ).

4.2 General case of the orthogonal FDSLRM

First of all, we formulate the basic assumptions required in our investigations on the structure and distribution of FDSLRM. As we mentioned in subsection 3, in real (practical) FDSLRM analysis we mainly work with orthogonal FDSLRMs (3.9). Therefore, in the rest of the paper we will focus primarily on this type of FDSLRMs.

The second assumption deals with a distribution of the FDSLRM observation $\boldsymbol{X}=(X(1),\ldots,X(n))^{\prime},$ $n\in\mathds{N}$. For now, we assume $\boldsymbol{X}$ satisfying (2.1) and having any absolutely continuous probability distribution with respect to some $\sigma$-finite measure.

Under these two assumptions, we can apply in the second stage of EBLUP-NE three estimation methods: natural estimators (NE), double ordinary least squares estimators (DOOLSE) and modified double ordinary least squares estimators (MDOOLSE).

Original natural estimators – NE

In the case of NE, we know the analytic solution of the estimation problem hancova_natural_2008 , which gives us always required non-negative estimates of FDSLRM variances. Employing results of hancova_natural_2008 , orthogonality condition (3.9) leads to the following form of NE

[TABLE]

where $\boldsymbol{e}=\boldsymbol{x}-\mathrm{F}\boldsymbol{\beta^{*}}=\mathrm{M_{F}}\boldsymbol{x}$ is nothing else than the vector of ordinary least squares (OLS) residuals in FDSLRM, $\mathrm{M_{V}}=\mathrm{I}_{n}-\mathrm{V}(\mathrm{V}^{\prime}\mathrm{V})^{-1}\mathrm{V}^{\prime}$ is the orthogonal projector onto the orthogonal element of the column space of $\mathrm{V}$ and $\boldsymbol{x}$ is an arbitrary realization of the FDSLRM observation $\boldsymbol{X}$. Vectors $\boldsymbol{v}_{j},j=1,2,\ldots,l$ are columns of design matrix $\mathrm{V}$. Moreover, in orthogonal FDSLRM BLUE $\boldsymbol{\beta}_{\boldsymbol{\upnu}}^{*}=(\mathrm{F}^{\prime}\mathrm{\Sigma_{\boldsymbol{\upnu}}^{-1}}\mathrm{F})^{-1}\mathrm{F}^{\prime}\mathrm{\Sigma_{\boldsymbol{\upnu}}^{-1}}\boldsymbol{x}$ from MME (3.2) is identical with ordinary least squares estimate $\boldsymbol{\beta^{*}}=(\mathrm{F}^{\prime}\mathrm{F})^{-1}\mathrm{F}^{\prime}\boldsymbol{x}$ not depending on $\boldsymbol{\upnu}$ stulajter_estimation_2004 .

As for computational complexity of NE, using elementary theory of complexity boyd_introduction_2018 , we get the complexity $2kn$ operations for the residual vector $\boldsymbol{e}$ and $4ln$ for NE from (4.1). Therefore, computing NE represents an algorithm with the complexity having order $\mathcal{O}(n)$ with respect to the realization length $n$.

Remark 4

For the purposes of software cross-checking of results and evaluating numerical precision of convex optimization algorithms, it is worth to formulate the calculation of NE as a convex optimization problem. Since NE employ the least-squares method, it is not complicated to show that in orthogonal FDSLRMs, NE can be obtained as a unique, always existing, non-negative solution of the following convex optimization problem

[TABLE]

where matrices in the objective function $f_{0}(\boldsymbol{\upnu})$ are defined by expressions

[TABLE]

In convex optimization, this kind of optimization problems belongs to the convex quadratic optimization problems or the norm approximation problems (boyd_convex_2009, , sec. 4.4, sec. 6.1).

Double least squares estimators – DOOLSE, MDOOLSE

From the optimization viewpoint, assuming general (not necessarily orthogonal) FDSLRM observation $\boldsymbol{X}$ (2.1), DOOLSE and MDOOLSE (stulajter_predictions_2002 , stulajter_estimation_2004 ) can be viewed as the following optimization problems for $\boldsymbol{\upnu}$ at given OLS residuals $\boldsymbol{e}$ (agrawal_rewriting_2018 , boyd_convex_2009 )

[TABLE]

where $\boldsymbol{e}$ are again OLS residuals as in the case of NE and $\mathrm{\Sigma_{\boldsymbol{\upnu}}}=\sigma_{0}^{2}\mathrm{I}_{n}+\mathrm{V}\mathrm{D}\mathrm{V}^{\prime}$ is the covariance matrix of $\boldsymbol{X}$. We call matrix $\mathrm{S}_{\boldsymbol{e}}\equiv\boldsymbol{e}\boldsymbol{e}^{\prime}$ matrix of residual products or more compactly the residual products matrix.

Remark 5

In the next text, some theoretical results can be written in one form for both of DOOLSE and MDOOLSE problems. To emphasize this fact, we will use one abbreviation with parenthesis $(\leavevmode\nobreak\ )$, e.g (M)DOOLSE, assuming that given results or considerations hold for both problems.

Using the geometrical language of Hilbert spaces (brockwell_time_2009 , stulajter_predictions_2002 ), (M)DOOLSE determine a covariance matrix with the parametric structure as in $\mathrm{\Sigma_{\boldsymbol{\upnu}}}$ but having the smallest Euclidean distance to the matrix $\mathrm{S}_{\boldsymbol{e}}$. In such case, (M)DOOLSE for orthogonal FDSLRMs could be computed geometrically as the orthogonal projection of $\mathrm{S}_{\boldsymbol{e}}$ onto the linear span $\mathscr{L}(\mathrm{V}_{0},\mathrm{V}_{1},\ldots,\mathrm{V}_{l})$ using the Gram matrix $\mathrm{G}$ of $\{\mathrm{V}_{0},\mathrm{V}_{1},\ldots,\mathrm{V}_{l}\}$ generated by the inner product $(\bullet,\bullet)=\mathrm{tr}(\bullet\cdot\bullet)$. Matrices $\mathrm{V}_{j},j=1,\ldots,l$ are given by (4.6), for DOOLSE: $\mathrm{V}_{0}=\mathrm{I}_{n}$ and for MDOOLSE: $\mathrm{V}_{0}=\mathrm{M_{F}}$.

According to stulajter_estimation_2004 , gajdos_kriging_2017 , the mentioned projection is given by

[TABLE]

where $\quad\mathrm{G}=\begin{pmatrix}n^{*}&\left\|\boldsymbol{v}_{1}\right\|^{2}&\left\|\boldsymbol{v}_{2}\right\|^{2}&\ldots&\left\|\boldsymbol{v}_{l}\right\|^{2}\\ \left\|\boldsymbol{v}_{1}\right\|^{2}&\left\|\boldsymbol{v}_{1}\right\|^{4}&0&\ldots&0\\ \left\|\boldsymbol{v}_{2}\right\|^{2}&0&\left\|\boldsymbol{v}_{2}\right\|^{4}&\ldots&0\\ \vdots&\vdots&\vdots&\ldots&\vdots\\ \left\|\boldsymbol{v}_{l}\right\|^{2}&0&0&\ldots&\left\|\boldsymbol{v}_{l}\right\|^{4}\end{pmatrix}\succ 0$,$\qquad\boldsymbol{q}=\begin{pmatrix}\boldsymbol{e}^{\prime}\boldsymbol{e}\\ (\boldsymbol{e}^{\prime}\boldsymbol{v}_{1})^{2}\\ (\boldsymbol{e}^{\prime}\boldsymbol{v}_{2})^{2}\\ \vdots\\ (\boldsymbol{e}^{\prime}\boldsymbol{v}_{l})^{2}\end{pmatrix},\\$ where for DOOLSE: $n^{*}\equiv n$, for MDOOLSE: $n^{*}\equiv n-k$, $\boldsymbol{e}$ are again OLS residuals, $\boldsymbol{v}_{j}$ are columns of $\mathrm{V}$.

However, the projection method can produce estimates out of the parametric space $\Upupsilon$ at a relatively high probability gajdos_kriging_2017 , i.e. in many cases it produces negative estimates. The proposed method (4.15) is not able to handle constraints given by $\Upupsilon$ reliably and in some simple way.

To distinguish the projection estimates from real DOOLSE (4.10), MDOOLSE (4.14) which are always non-negative, we suggest to be more accurate and use full names for real DOOLSE and MDOOLSE: non-negative DOOLSE (NN-DOOLSE) and non-negative MDOOLSE (NN-MDOOLSE) This specific way is similarly used e.g. in the case of Rao’s MINQUE (rao_estimation_1988, , chap. 5). For the projection-based DOOLSE and MDOOLSE, without considering nonnegativity constraints, we leave the original acronyms DOOLSE, MDOOLSE.

Applying the basic convex quadratic optimization theory (cornuejols_optimization_2018 , bertsekas_convex_2009 ) in orthogonal FDSLRM, we can rewrite NN-(M)DOOLSE as strictly convex quadratic problems, whose solutions $\tilde{\boldsymbol{\upnu}}$ always exist and are unique global minimizers on $\Upupsilon^{*}$ (see Propositions 2 and 4 in the Appendix)

[TABLE]

The most important result of convex optimization for NN-(M)DOOLSE problems are fundamental Karush-Kuhn-Tucker (KKT) optimality conditions, which in orthogonal FDSLRM provide a necessary and sufficient condition for optimality of solutions (see Proposition 3 in the Appendix). In practice, finding solutions of convex optimization problems is nothing else than solving the corresponding KKT conditions analytically or in general numerically.

In our case of NN-(M)DOOLSE, the KKT conditions can be represented by a set of $2^{l}$ nonsingular linear systems of equations, where each of them is given by matrix $\mathrm{K}$ derived from $\mathrm{G}$ and by vector $\boldsymbol{q}$ from (4.15). Our theoretical results in the Appendix show that only one of the linear systems always gives us all required non-negative elements of $\tilde{\boldsymbol{\upnu}}$.

On this basis, we can establish a simple KKT optimization algorithm which run through given linear systems, compute their analytical solution and stops when the solution appears non-negative. In other words, the proposed KKT algorithm always founds the required NN-(M)DOOLSE $\tilde{\boldsymbol{\upnu}}_{+}$ in at most $2^{l}$ steps. The algorithm is summarized in tab. 3, whereas the proof and details are in the Appendix.

As for computational complexity of the KKT optimization algorithm, all calculations of input have complexity $\mathcal{O}(n)$, whereas the body of the algorithm (steps 1-3) is $\mathcal{O}(l^{2}\cdot 2^{l})$. Since fixed number $l$ of variance parameters in matrix $\mathrm{D}_{\boldsymbol{\upnu}}={Cov}\{\boldsymbol{Y}\}$ is usually much smaller than $n$ (typically $l^{2}\cdot 2^{l}$ is of $\mathcal{O}(n)$), then the complete algorithm has the leading order $\mathcal{O}(n)$ with respect to $n$.

It is the same as in the case of NE and generally typical only for computationally fastest analytical solutions of the KKT conditions.

Computational tools for orthogonal FDSLRM

In nonlinear optimization, there is a variety of highly efficient, fast and reliable open-source and commercial software packages (cornuejols_optimization_2018, , chap. 20). We have chosen one of the most well-known open source libraries for solving convex optimization tasks – CVXPY diamond_cvxpy:_2016 and its R version CVXR fu_cvxr:_2019 .

CVXPY is a scientific Python library but also a language with very simple, readable syntax not requiring any expertise in convex optimization and its PC implementation. CVXPY allows the user to specify the mathematical optimization problem naturally following normal mathematical notation as we can see in computing NN-DOOLSE (fig. 1) where a code easily mimics the non-negative DOOLSE mathematical formulation (4.10) with $\mathrm{\Sigma_{\boldsymbol{\upnu}}}=\sigma_{0}^{2}\mathrm{I}_{n}+\mathrm{V}\mathrm{D}\mathrm{V}^{\prime}=\upnu_{0}\mathrm{I}_{n}+\mathrm{V}\mathrm{diag}\left\{\upnu_{j}\right\}\mathrm{V}^{\prime}$.

CVXPY, implements not only convex optimization solvers using interior-point numerical methods which are extremely reliable and fast, but also verifies convexity of the given problem using rules of disciplined convex programming grant_disciplined_2006 . In FDSLRM, interior-points numerical methods have complexity $n^{3}$ for one iteration or $\log(1/\varepsilon)$ times bigger for the complete computation with a precision $\varepsilon$ of the required optimal solution boyd_convex_2009 .

Remark 6

It is worth to mention that both CVXPY and CVXR are able to solve non-negative DOOLSE and MDOOLSE in a general FDSLRM without the orthogonality condition. CVXPY was inspired by MATLAB optimization package CVX grant_cvx:_2018 still used in many references, e.g. cornuejols_optimization_2018 . However, in CVXPY (CVXR) the user can easily combine convex optimization and simplicity of Python (R) language together with its high-level features such as object-oriented design or parallelism.

4.3 Gaussian orthogonal FDSLRM

In the LMM framework, the most usual $\boldsymbol{X}$ distribution used in practice is represented by the multivariate normal (Gaussian) distribution. Therefore, in addition to the orthogonality of FDSLRM, we will require in this section the distributional assumption $\boldsymbol{X}\sim\mathcal{N}_{n}(\mathrm{F}\boldsymbol{\beta},\mathrm{\Sigma_{\boldsymbol{\upnu}}})$. We will refer to FDSLRM under these assumptions as to Gaussian orthogonal FDSLRM.

In the case of Gaussian orthogonal FDSLRM, we can add to our investigation last two estimation methods from tab. 2: maximum likelihood estimators (MLE) and residual maximum likelihood estimators (REMLE).

Maximum likelihood estimators – MLE, REMLE

Both MLE and REMLE of variances $\boldsymbol{\upnu}$ provide estimates maximizing ML and REML loglikelihood functions (logarithms of likelihoods). Using simple and clear arguments from Lamotte’s paper lamotte_direct_2007 , we can easily form the ML loglikelihood function $l_{\small\mathrm{m}}(\boldsymbol{\upnu},\boldsymbol{x})$ and REML loglikelihood function $l_{\small\mathrm{r}}(\boldsymbol{\upnu},\boldsymbol{x})$, assuming general (not necessarily orthogonal) FDSLRM observation $\boldsymbol{X}$ (2.1), as

[TABLE]

where $\boldsymbol{e}_{\boldsymbol{\upnu}}=\boldsymbol{x}-\mathrm{F}\boldsymbol{\beta}_{\boldsymbol{\upnu}}^{*}$ are FDSLRM residuals, known as marginal residuals in LMM singer_graphical_2017 , $\boldsymbol{x}$ is an arbitrary realization of the FDSLRM observation $\boldsymbol{X}$ with covariance matrix $\mathrm{\Sigma_{\boldsymbol{\upnu}}}$, $\boldsymbol{\beta}_{\boldsymbol{\upnu}}^{*}$ is BLUE of $\boldsymbol{\beta}$ in MME (3.2) and $\left\|\bullet\right\|_{\mathrm{\Sigma_{\boldsymbol{\upnu}}^{-1}}}^{2}=(\bullet)^{\prime}\mathrm{\Sigma_{\boldsymbol{\upnu}}^{-1}}(\bullet)$ is a generalized vector norm given by $\mathrm{\Sigma_{\boldsymbol{\upnu}}^{-1}}$.

Remark 7

In Lammote’s paper lamotte_direct_2007 , which presents a direct derivation of REML likelihood function in LMM using familiar linear algebra operations, we can find that the formulation of ML and REML likelihood in LMM under the assumption of normality was originally done in Harville’s works harville_bayesian_1974 , harville_maximum_1977 . However, LMM references up to 2007 show the explicit form of REML likelihood for LMM rarely. In Štulajter’s FDSLRM reference works dealing with MLE stulajter_predictions_2002 , stulajter_estimation_2004 REML likelihood for FDSLRM is also absent. According to Lammote lamotte_direct_2007 , the reason why seems to be Harville’s too sophisticated, difficult and indirect derivation.

Maximum likelihood and residual maximum likelihood estimation of variance components in LMMs produce, in general, no analytical expressions for the estimators (searle_variance_2009, , chap. 6). They exist only for $\upnu_{0}>0$. According to (demidenko_mixed_2013, , sec. 2.5, thm 4) MLE and REMLE estimates do not exist in general FDSLRM with LMM observation $\boldsymbol{X}$ (2.1) if and only if a realization $\boldsymbol{x}$ of $\boldsymbol{X}$ belongs to $\mathscr{L}(\mathrm{F},\mathrm{V})$, linear span of columns $\mathrm{F}$ and $\mathrm{V}$. This condition is also equivalent with $\boldsymbol{e}\in\mathscr{L}(\mathrm{V})$.

In orthogonal FDSLRM, the MLE and REMLE can be rewritten as the optimization problems for $\boldsymbol{\upnu}$ (agrawal_rewriting_2018 , boyd_convex_2009 ) at given OLS residuals $\boldsymbol{e}$ (for orthogonal FDSLRM $\boldsymbol{e}=\boldsymbol{e}_{\boldsymbol{\upnu}}$, stulajter_estimation_2004 )

[TABLE]

[TABLE]

where $\left\|\boldsymbol{e}\right\|_{\mathrm{\Sigma_{\boldsymbol{\upnu}}^{-1}}}^{2}=\boldsymbol{e}^{\prime}\mathrm{\Sigma_{\boldsymbol{\upnu}}^{-1}}\boldsymbol{e}=\mathrm{tr}\left(\mathrm{S}_{\boldsymbol{e}}\mathrm{\Sigma_{\boldsymbol{\upnu}}^{-1}}\right)$ and $\mathrm{S}_{\boldsymbol{e}}$ is the residual products matrix used in definitions (4.10), (4.14) of (M)DOOLSE.

Concerning convex optimization, (RE)MLE problems (4.29), (4.25) are generally not convex problems. However in orthogonal FDSLRM, according to (boyd_convex_2009, , sec. 4.1.3), using an appropriate bijective transformation, we can reformulate (RE)MLE as equivalent convex optimization problems (see Proposition 5).

Employing the theory of convex optimization, we can prove that equivalent (RE)MLE convex problems are strictly convex, for $\boldsymbol{e}\not\in\mathscr{L}(\mathrm{V})$ have always a unique solution, unique global minimizer satisfying corresponding KKT conditions (see Propositions 7, 6 and their proofs in the Appendix). These conditions are again necessary and sufficient for optimal solutions.

Combining all our theoretical results, we arrive to the main theoretical result of the section. It can be shown that for $\upnu_{0}\neq 0$ (or $\boldsymbol{e}\not\in\mathscr{L}(\mathrm{V})$) the KKT conditions for NN-(M)DOOLSE are equivalent with the KKT conditions for (RE)MLE and the bijective transformation. Or in other words, we can formulate the following theorem.

Theorem 4.1 (equivalence between NN-(M)DOOLSE and (RE)MLE)

In a Gaussian orthogonal FDSLRM non-negative (M)DOOLSE are almost sure equal to (RE)MLE, i.e.

[TABLE]

Proof

See the Appendix.

Computational tools for Gaussian orthogonal FDSLRM

The most important and useful consequence of Theorem 4.1 is that we can compute MLE or REMLE in orthogonal FDSLRM as quadratic NN-(M)DOOLSE. Therefore, we can use our KKT algorithm of order $\mathcal{O}(n)$ or computational tools CVXPY (CVXR) with $\mathcal{O}(n^{3})$ described in the previous section. Moreover, for $\boldsymbol{e}\in\mathscr{L}(\mathrm{V})$ NN-(M)DOOLSE do not fail as (RE)MLE, but naturally extend them.

Since a Gaussian orthogonal FDSLRM is also Gaussian LMM, we can also apply current LMM packages. In our previous paper gajdos_kriging_2017 , we stated that no current package in R is directly and effectively suitable for FDSLRM. Thanks to a detailed study of demidenko_mixed_2013 , galecki_linear_2013 we were successfully directed and instructed how to implement the FDSLRM variance structure into one of the best R packages for LMM known as nlme (pinheiro_mixed-effects_2009 , pinheiro_nlme:_2018 ). After that, inspired by nlme, we found another R package called sommer, also suitable for FDSLRM fitting (covarrubias-pazaran_genome-assisted_2016 , cornuejols_optimization_2018 ).

Simultaneously, thanks to online computing environment CoCalc for SageMath stein_sage_2019 with the possibility to run computations and codes for free in many other programming languages and open softwares, we are also able to run and test MATLAB function mixed (witkovsky_mixed_2018 , witkovsky_estimation_2012 ) primarily intended for estimating variances in LMMs. On the basis of successful tests, we wrote its R version called MMEinR, which is also available in our GitHub Repository gajdos_mmeinr:_2019 .

Remark 8

The mentioned LMM packages can also handle (RE)MLE in the general FDSLRM without orthogonality restriction. The packages use iterative methods based on EM algorithm or Henderson MME together with some version of the Newton-Raphson method whose complexity is generally at least $n^{3}$.

We are fully aware that another efficient implementation for estimating variances in LMMs provides lme4 package bates_fitting_2015 or SAS package PROC MIXED stroup_sas_2018 , sasinstituteinc._sas/stat_2018 . As for lme4, which is much faster than nlme, we found that the package does not allow implementing FDSLRM using lme4 standard input procedure.

Finally, as SAS laymans, we also tried a university edition of SAS, free from 2014, but we left it after running into initial problems with Windows installation process in a virtual machine environment and subsequently with a more sophisticated programming language than Python or R. Therefore our knowledge with SAS remained only at the theoretical level.

5 Application in real data examples

5.1 Three real data sets and their FDSLRMs

We illustrate the obtained theoretical results and the performance of the proposed EBLUP-NE method of variances $\boldsymbol{\upnu}$ in the FDSLRM (1.1) using three real data sets: electricity consumption, tourism and cyber security. For all data sets, we identified the most parsimonious structure of the FDSLRM using an iterative process of the model building and selection based on exploratory tools of spectral analysis and LMM theory. Details of our analysis and modeling can be found in our easy reproducible Jupyter notebooks freely available at our GitHub repository gajdos_fdslrm_2019 .

Electricity consumption

As the first real data example, we have the econometric time series data set, representing hourly observations of the consumption of the electric energy (in kWh) in an department store. The number of time series observations is $n=24$. The data was adapted from stulajter_estimation_2004 . The consumption data can be fitted by the following Gaussian orthogonal FDSLRM stulajter_estimation_2004 :

[TABLE]

with $k=3,l=4,\;\boldsymbol{Y}=(Y_{1},Y_{2},Y_{3},Y_{4})^{\prime}\sim\mathcal{N}_{4}(\boldsymbol{0}_{4},\mathrm{D}_{\boldsymbol{\upnu}}),\;w(t)\sim iid\,\mathcal{N}(0,\sigma_{0}^{2})$ and $\boldsymbol{\upnu}=(\sigma_{0}^{2},\sigma_{1}^{2},\sigma_{2}^{2},\sigma_{3}^{2},\sigma_{4}^{2})^{\prime}\in\mathds{R}^{5}_{+}$. Frequencies $\omega_{1},\omega_{2},\omega_{3}$ are suitable Fourier frequencies from the time series periodogram in spectral analysis (priestley_spectral_2004 , stulajter_estimation_2004 , gajdos_kriging_2017 ).

Since the time series data set contains only 24 observations and we have to estimate three regression and five variance parameters of the FDSLRM, this real data FDSLRM example should be considered only as a toy example. Taking two sets of Fourier frequencies $\left(2\pi/24,2\pi\cdot 3/24,2\pi\cdot 4/24\right)$ or $\left(2\pi/24,2\pi\cdot 2/24,2\pi\cdot 3/24\right)$, we get two toy models previously introduced in gajdos_kriging_2017 and stulajter_estimation_2004 . These models give us the opportunity to check our numerical results and to demonstrate how in principle FDSLRM estimation methods and computational tools work.

Tourism

In this econometric FDSLRM application, we consider the time series data set, called visnights, representing total quarterly visitor nights (in millions) from 1998-2016 in one of the regions of Australia – inner zone of Victoria state. The number of time series observations is $n=76$. The data was adapted from hyndman_forecasting:_2018 .

The Gaussian orthogonal FDSLRM fitting the tourism data has the following form (for details see our Jupyter notebook tourism.ipynb):

[TABLE]

with $k=3,l=3,\;\boldsymbol{Y}=(Y_{1},Y_{2},Y_{3})^{\prime}\sim\mathcal{N}_{3}(\boldsymbol{0}_{3},\mathrm{D}_{\boldsymbol{\upnu}})\,$ and $\boldsymbol{\upnu}=(\sigma_{0}^{2},\sigma_{1}^{2},\sigma_{2}^{2},\sigma_{3}^{2})^{\prime}\in\mathds{R}^{4}_{+}$.

Cyber attacks

Our final FDSLRM application describes the real time series data set representing the total weekly number of cyber attacks against a honeynet – an unconventional tool which mimics real systems connected to Internet, like business or school computers intranets, to study methods, tools and goals of cyber attackers. Data, taken from sokol_prediction_2018 , were collected from November 2014 to May 2016 in CZ.NIC honeynet consisting of Kippo honeypots in medium-interaction mode. The number of time series observations is $n=72$.

The suitable FDSLRM, after a preliminary logarithmic transformation of data $Z(t)=\log X(t)$, is again Gaussian orthogonal (for details see our Jupyter notebook cyberattacks.ipynb) and in comparison with previous models (5.30), (5.31) has the simplest structure:

[TABLE]

with $k=3,l=2,\;\boldsymbol{Y}=(Y_{1},Y_{2})^{\prime}\sim\mathcal{N}_{2}(\boldsymbol{0}_{2},\mathrm{D}_{\boldsymbol{\upnu}})\,$ and $\boldsymbol{\upnu}=(\sigma_{0}^{2},\sigma_{1}^{2},\sigma_{2}^{2})^{\prime}\in\mathds{R}^{3}_{+}$.

5.2 Numerical results of estimating $\boldsymbol{\upnu}$

For cross-checking purposes, we realized our numerical computations in both Python and R based software tools (or packages). As we described in previous sections, we implemented own algorithms and methods in SciPy, SageMath and R, particularly analytical expressions (4.1) for NE and the KKT algorithm for NN-(M)DOOLSE (tab. 3) in orthogonal FDSLRM.

Simultaneously, we confirmed the same estimations using CVXPY (or CVXR) package based on convex optimization and analogically using up-to-date standard LMM R packages nlme, MMEinR and sommer.

Detailed computational results for all three data sets using all corresponding relevant tools can be found and reproduced in our collection of Jupyter notebooks (asterisk * represents a name of data and specific computational tool): PY-estimation-.ipynb and R-estimation-.ipynb.

Table 4 summarizes all types of considered initial estimates and their corresponding EBLUP-NE in all four models. The computations also confirmed Theorem 4.1 therefore we wrote results for MLE and REMLE only.

Finally, we point out that thanks to SageMath and our very fast KKT algorithm, we were able to compute (in real time) NN-(M)-DOOLSE in toy examples with infinite precision – as the exact closed-form (algebraic) numbers. To get an explicit idea, e.g. for our first toy model, we got these results for NN-MDOOLSE $\tilde{\boldsymbol{\upnu}}_{+}$ identical with REMLE

[TABLE]

It means that we can compute real errors in results for all used computational tools to explore their quality from the viewpoint of numerical precision. At the same time, we also watched a run time for particular computational tools using Jupyter notebook extension ExecuteTime, which provided a preliminary comparison of real execution times.

6 Conclusions

We suggested and investigated an alternative, new method EBLUP-NE based on empirical BLUPs for estimating variances in time series modeled by FDSLRM. The estimation method can be also viewed, analogously like EBLUP (rao_small_2015 , witkovsky_estimation_2012 ), as a two-stage iterative method with one step in the iteration.

EBLUP-NE are invariant quadratic, non-negative estimators whose simple computational form and subsequent first and second-moment statistical properties are given by two special matrices – the Schur complement and Gram matrix zhang_schur_2005 . The method can be used not only in the case of normally distributed time series data, but for any absolutely continuous probability distribution of time series data.

As initial starting estimates for EBLUP-NE, we can principally employ any of the previously used methods based on least squares (NE, NN-DOOLSE, NN-MDOOLSE) or maximum likelihood (MLE, REMLE). Applications of FDSLRM with the EBLUP-NE on three real data sets (electricity consumption, tourism, cyber attacks) providing two toy and two real models indicate that the method computationally gives at least comparable results with REMLE or NN-MDOOLSE, but in faster run time (approximately 10-1000 times on the standard PC).

Due to lack of computational implementation for FDSLRM modeling, which would be generally available and readily applicable, and the fact that FDSLRM least squares and maximum likelihoods estimation procedures are more than 10 years old, we revisited and updated theoretical and computational knowledge dealing with these methods. Specifically, applying the convex optimization theory, we reformulated all estimation methods as convex optimization problems in the so-called orthogonal FDSLRM, the most usual form of FDSLRM used in practice.

We formulated the KKT optimality conditions, which, unlike likelihood equations, are necessary and sufficient conditions for optimal solutions of (RE)MLE or NN-(M)DOOLSE on extended parametric space $\Upupsilon=[0,\infty)^{l+1}$. KKT optimality conditions dictate not only the exact existence conditions of estimates but they also solve the well-known problem dealing with standardly used likelihood equations (christensen_plane_2011, , chap. 12), (searle_variance_2009, , chap. 6) in LMM where their solutions for (RE)MLE or NN-(M)DOOLSE may not be required estimates — they may be out of the parameter space or they may be other than the maximum or minimum.

Moreover, using KKT optimality conditions in orthogonal FDSLRM, we proved the equivalence of NN-(M)DOOLSE and (RE)MLE with propability 1. This most important theoretical result of the paper is the stronger and more general result than in stulajter_estimation_2004 proved only for interior points of $\Upupsilon$.

Simultaneously, the convex optimization theory brought us to the new KKT algorithm for computing NN-(M)DOOLSE, equivalent to (RE)MLE, with double floating-point precision $\epsilon<10^{-15}$ as the default precision of outputs and with computational complexity $\mathcal{O}(n)$. Such an algorithm which we implemented in SciPy, SageMath and R, which at the default precision level is $10^{7}$ times more accurate and approximately $n^{2}$ faster than the best current Python- or R-based standard computational tools, can be used in effective computational time series research to study properties of FDSLRM (Monte Carlo and bootstrap methods, kreiss_bootstrap_2012 ).

Regarding computational aspects, we were also successful in the identification and demonstration of consistent results for the real data applications of FDSLRM in several free, open-source current standard computational tools — namely CVXPY, CVXR (R version of CVXPY) packages for convex optimization and LMM R packages nlme, sommer and MMEinR (our R version of Witkovsky’s MATLAB mixed function). These results and procedures can be freely viewed in our 16 Jupyter notebooks which are easily readable, sharable, reproducible and modifiable directly in our GitHub repository (gajdos_fdslrm_2019 ). Open-source Jupyter technology with Python and R packages also solved our problem stated in gajdos_kriging_2017 that no current package in R was directly and effectively suitable for FDSLRM.

Finally, our investigation has also brought new questions for further research. There is definitely a need for more exact and detailed analysis based on a simulation study and general EBLUP theory (zadlo_mse_2009 , witkovsky_estimation_2012 , rao_small_2015 ) focusing on EBLUP-NE quality with respect to previously used estimation methods and the performance of the EBLUP-NE method itself with different initial starting points, using different computational tools in various probability distributions. These research questions are under our investigation and will be published in the near future.

In connection with our current computational research in FDSLRM, but also for real time series data analysis and forecasting, we started to build our own R package (see a preliminary, fully functional version at gajdos_fdslrmR_2019 ) on mentioned LMM R packages to manipulate readily with FDSLRM concepts and procedures.

Finally, our results in the paper can be seen reciprocally as contributions to convex optimization and LMM methodology. Particularly, our convex optimization application in the context of time series modeling has become another one from a wide variety of application areas of convex optimization boyd_convex_2009 . Since FDSLRM describing $n$ observed time series values is also a special type of LMM, our EBLUP-NE and the very fast, accurate KKT algorithm to compute (RE)MLE may have potential to be used in computational research and applications dealing with LMMs.

Acknowledgements.

We are very grateful especially to František Štulajter (Comenius University in Bratislava, SK) for his guidance during our early contact and study of FDSLRM concepts and to Viktor Witkovský (Slovak Academy of Science, SK) for his recommendations and deep insights dealing with the LMM methodology and its computational tools connected to FDSLRM. We would also like to thank Ivan Žezula and Daniel Klein (both from P. J. Šafárik University in Košice, SK) for a valuable feedback dealing with considered estimation methods. Concerning applied computational tools, we would also like to acknowledge involvement and recommendations of Giovanny Covarrubias-Pazaran (University of Wisconsin, US) in using sommer package, Steve Diamond and Stephen Boyd (Stanford University, US) in using CVXPY and CVXR. Finally, in connection with using SageMath, Jupyter and GitHub, our specials thanks go namely to William Stein (University of Washington, US), Harald Schilly (Universität Wien, AT), Erik Bray (Université Paris-Sud, FR), Luca de Feo (Université de Versailles Saint-Quentin-en-Yvelines, FR) and Charles J. Weiss (Augustana University, US).

Appendix

Acronyms and abbreviations

The BLUP-NE method

Proof (Theorem 3.1)

Employing Lemma 1 and the well-known standard expressions for mean values and covariances of invariant quadratic estimators (see e.g. christensen_plane_2011 , puntanen_formulas_2013 )) $E_{\boldsymbol{\upnu}}\left\{\boldsymbol{X}^{\prime}\mathrm{A}\boldsymbol{X}\right\}=\mathrm{tr}(\mathrm{A}\mathrm{\Sigma_{\boldsymbol{\upnu}}})$ and if $\boldsymbol{X}\sim\mathcal{N}_{n}(\mathrm{F}\boldsymbol{\beta},\mathrm{\Sigma_{\boldsymbol{\upnu}}})$ then ${Cov}_{\boldsymbol{\upnu}}\left\{\boldsymbol{X}^{\prime}\mathrm{A}\boldsymbol{X},\boldsymbol{X}^{\prime}\mathrm{B}\boldsymbol{X}\right\}=2\,\mathrm{tr}(\mathrm{A}\mathrm{\Sigma_{\boldsymbol{\upnu}}}\mathrm{B}\mathrm{\Sigma_{\boldsymbol{\upnu}}})$, we have

1. $(0)$

$E_{\boldsymbol{\upnu}}\left\{\mathring{{\sigma}}_{j}^{2}\right\}=E(\boldsymbol{X}\mathbb{t}_{j}\mathbb{t}_{j}^{\prime}\boldsymbol{X})=\mathrm{tr}(\mathbb{t}_{j}\mathbb{t}_{j}^{\prime}\mathrm{\Sigma_{\boldsymbol{\upnu}}})=\mathrm{tr}(\mathbb{t}_{j}^{\prime}\mathrm{\Sigma_{\boldsymbol{\upnu}}}\mathbb{t}_{j})=(\mathbb{T}\mathrm{\Sigma_{\boldsymbol{\upnu}}}\mathbb{T}^{\prime})_{jj}$.

According to (3) of Lema 1 $(\mathbb{T}\mathrm{\Sigma_{\boldsymbol{\upnu}}}\mathbb{T}^{\prime})_{jj}=\mathrm{D}_{jj}-\sigma^{2}(\mathbb{W}^{-1})_{jj}=\sigma^{2}_{j}-\sigma^{2}(\mathbb{W}^{-1})_{jj}$. 2. $(0)$

$D_{\boldsymbol{\upnu}}\{\mathring{{\sigma}}_{j}^{2}\}=2\mathrm{tr}(\mathbb{t}_{j}\mathbb{t}_{j}^{\prime}\mathrm{\Sigma_{\boldsymbol{\upnu}}}\mathbb{t}_{j}\mathbb{t}_{j}^{\prime}\mathrm{\Sigma_{\boldsymbol{\upnu}}})=2\mathrm{tr}(\mathbb{t}_{j}^{\prime}\mathrm{\Sigma_{\boldsymbol{\upnu}}}\mathbb{t}_{j}^{\prime}\mathbb{t}_{j}\mathrm{\Sigma_{\boldsymbol{\upnu}}}\mathbb{t}_{j})=2(\mathbb{T}\mathrm{\Sigma_{\boldsymbol{\upnu}}}\mathbb{T}^{\prime})_{jj}$

As a consequence of (1) $D_{\boldsymbol{\upnu}}\{\mathring{{\sigma}}_{j}^{2}\}=2(\sigma^{2}_{j}-\sigma^{2}(\mathbb{W}^{-1})_{jj})^{2}$. 3. $(0)$

In a similar way as in (2)

${Cov}_{\boldsymbol{\upnu}}\left\{\mathring{{\sigma}}_{i}^{2},\mathring{{\sigma}}_{j}^{2}\right\}=2\mathrm{tr}(\mathbb{t}_{i}\mathbb{t}_{i}^{\prime}\mathrm{\Sigma_{\boldsymbol{\upnu}}}\mathbb{t}_{j}\mathbb{t}_{j}^{\prime}\mathrm{\Sigma_{\boldsymbol{\upnu}}})=2\mathrm{tr}(\mathbb{t}_{i}^{\prime}\mathrm{\Sigma_{\boldsymbol{\upnu}}}\mathbb{t}_{j}\mathbb{t}_{j}^{\prime}\mathrm{\Sigma_{\boldsymbol{\upnu}}}\mathbb{t}_{i})=2(\mathbb{t}_{i}^{\prime}\mathrm{\Sigma_{\boldsymbol{\upnu}}}\mathbb{t}_{j})^{2}=\\ =2(\mathbb{T}\mathrm{\Sigma_{\boldsymbol{\upnu}}}\mathbb{T}^{\prime})^{2}_{ij}=2(0-\sigma^{2}(\mathbb{W}^{-1})_{ij})^{2}=2(\sigma^{2}(\mathbb{W}^{-1})_{ij})^{2}$. 4. $(0)$

${MSE}_{\boldsymbol{\upnu}}\{\mathring{{\sigma}}_{j}^{2}\}=E_{\boldsymbol{\upnu}}\left\{{(\mathring{{\sigma}}_{j}^{2}-\sigma^{2}_{j})}^{2}\right\}={(E_{\boldsymbol{\upnu}}\left\{\mathring{{\sigma}}_{j}^{2}\right\}-\sigma^{2}_{j})}^{2}+D_{\boldsymbol{\upnu}}\{\mathring{{\sigma}}_{j}^{2}\}=\\ =(\sigma^{2}_{j}-\sigma^{2}(\mathbb{W}^{-1})_{jj}-\sigma^{2}_{j})^{2}+2(\sigma^{2}_{j}-\sigma^{2}(\mathbb{W}^{-1})_{jj})^{2}=\\ =\sigma^{4}(\mathbb{W}^{-1})^{2}_{jj}+2(\sigma^{2}_{j}-\sigma^{2}(\mathbb{W}^{-1})_{jj})^{2}$.

∎

Existence of inversions for $\mathbb{U}$, $\mathbb{H}_{\mathrm{V}}$

To prove nonsingularity of $\mathbb{U}$, $\mathbb{H}_{\mathrm{V}}$, it is sufficient to show that both matrices have non-zero determinants. Using idempotence of $\mathrm{M_{F}}$ and expression (searle_matrix_2017, , sec. 6.8)

[TABLE]

we can write for determinants of $\mathbb{U}$, $\mathbb{H}_{\mathrm{V}}$

[TABLE]

Now we can see that the sum of positive definite matrix $\sigma_{0}^{2}\mathrm{I}_{l}\,\,(\sigma_{0}>0)$ and each of positive semidefinite matrices $\mathrm{V}\mathrm{D}\mathrm{V}^{\prime},\,\,\mathrm{M_{F}}\mathrm{V}\mathrm{D}\mathrm{V}^{\prime}\mathrm{M_{F}}$ is a positive definite matrix, whose determinant is always positive.

Double least squares estimators NN-(M)DOOLSE

Using the standard inner product of matrices defined as $(\bullet,\bullet)=\mathrm{tr}(\bullet,\bullet)$, generating the Euclidean norm $\left\|\bullet\right\|$ of a matrix, and basic properties of the trace function, we can easily rewrite expressions (4.10), (4.14) for the objective functions as quadratic forms. The particular form of these quadratic forms in optimization problems (M)DOOLSE is described by the following proposition.

Proposition 2 (NN-(M)DOOLSE as quadratic optimization problems)

In an orthogonal FDSLRM the NN-(M)DOOLSE problems are convex quadratic optimization problems in the form

[TABLE]

\qquad\boldsymbol{q}=\left(\begin{array}[]{c}\boldsymbol{e}^{\prime}\boldsymbol{e}\\ (\boldsymbol{e}^{\prime}\boldsymbol{v}_{1})^{2}\\ (\boldsymbol{e}^{\prime}\boldsymbol{v}_{2})^{2}\\ \vdots\\ (\boldsymbol{e}^{\prime}\boldsymbol{v}_{l})^{2}\end{array}\right)*, \quad\mathrm{G}=\left(\begin{array}[]{ccccc}n^{*}&\left\|\boldsymbol{v}_{1}\right\|^{2}&\left\|\boldsymbol{v}_{2}\right\|^{2}&\ldots&\left\|\boldsymbol{v}_{l}\right\|^{2}\\ \left\|\boldsymbol{v}_{1}\right\|^{2}&\left\|\boldsymbol{v}_{1}\right\|^{4}&0&\ldots&0\\ \left\|\boldsymbol{v}_{2}\right\|^{2}&0&\left\|\boldsymbol{v}_{2}\right\|^{4}&\ldots&0\\ \vdots&\vdots&\vdots&\ldots&\vdots\\ \left\|\boldsymbol{v}_{l}\right\|^{2}&0&0&\ldots&\left\|\boldsymbol{v}_{l}\right\|^{4}\end{array}\right)\succ 0

where for *NN-DOOLSE: * $n^{*}\equiv n$,  *NN-MDOOLSE: $n^{*}\equiv n-k$,

According to the fundamental KKT theorem of convex optimization (boyd_convex_2009, , chap.5) which handles convex optimization problems with constraints in the form of inequalities giving a list of the so-called Karush-Kuhn-Tucker (KKT) optimality conditions, we consider our optimization problem NN-(M)DOOLSE in the Lagrange multiplier form with Lagrangian $L$ as a sum of the objective function $f_{0}(\boldsymbol{\upnu})$ and a linear combination of multipliers $\boldsymbol{\uplambda}=(\uplambda_{0},\uplambda_{1},\ldots,\uplambda_{l})^{\prime}\succeq\boldsymbol{0}_{l+1}$ and constraints $-\boldsymbol{\upnu}=(\upnu_{0},\upnu_{1},\ldots,\upnu_{l})^{\prime}\preceq\boldsymbol{0}_{l+1}$.

Then KKT optimality conditions for (6.37) become necessary and sufficient conditions of optimality, satisfied at any local optimal solution, being represented by three sets of conditions: (1) primal and dual feasibility: $\boldsymbol{\upnu}\succeq\boldsymbol{0}_{l+1}$, $\lambda\succeq\boldsymbol{0}_{l+1}$, (2) stationarity of the Lagrangian: $\nabla_{\boldsymbol{\upnu}}L(\boldsymbol{\upnu},\boldsymbol{\uplambda})=\boldsymbol{0}_{l+1}$ and (3) complementary slackness: $\upnu_{j}\uplambda_{j}=0,\,j\in\{0,1,\ldots,l\}$. Computing the gradient of the Lagrangian, we arrive at the following proposition.

Proposition 3 (KKT optimality conditions for NN-(M)DOOLSE)

In an orthogonal FDSLRM consider the non-negative NN-(M)DOOLSE problems in the Lagrange multiplier form with the Lagrangian

[TABLE]

Then, a necessary and sufficient condition for $\tilde{\boldsymbol{\upnu}},\tilde{\boldsymbol{\uplambda}}$ to be problems’ optimal solution is

1. $(0)$

$\boldsymbol{\upnu}\succeq\boldsymbol{0}_{l+1},\quad\boldsymbol{\uplambda}\succeq\boldsymbol{0}_{l+1}$** 2. $(0)$

$\begin{pmatrix}\mathrm{G}&-\mathrm{I}_{l+1}\\ \end{pmatrix}\begin{pmatrix}\boldsymbol{\upnu}\\ \boldsymbol{\uplambda}\end{pmatrix}=\boldsymbol{q}$** 3. $(0)$

$\boldsymbol{\upnu}\circ\boldsymbol{\uplambda}=\boldsymbol{0}_{l+1}\quad\left(\Leftrightarrow\upnu_{j}\uplambda_{j}=0,\,\,j\in\{0,1,\ldots,l\}\right)$.

As for the existence of optimal solutions, we can apply the well-known basic results for minimizing quadratic forms (bertsekas_convex_2009, , chap. 3). Since the Hessian of (6.37) equal to $\nabla_{\boldsymbol{\upnu}}^{2}f_{0}(\boldsymbol{\upnu})=2\mathrm{G}$ is a positive definite matrix, like $\mathrm{G}$, then $f_{0}(\boldsymbol{\upnu})$ is strictly convex and also coercive proper function. These two properties are sufficient conditions for NN-(M)DOOLSE (bertsekas_convex_2009, , Weierstrass’ theorem, p. 119) to have always a unique global optimal solution $\tilde{\boldsymbol{\upnu}}$.

Employing the familiar Bessel’s inequality $\sum\limits_{j=1}^{l}(\boldsymbol{e}^{\prime}\boldsymbol{v}_{j}/\left\|\boldsymbol{v}_{j}\right\|)^{2}\leq\left\|\boldsymbol{e}\right\|^{2}$, where equality holds if and only if $\boldsymbol{e}\in\mathscr{L}(\mathrm{V})$, it is also not difficult to see that KKT conditions (1-3) rewritten as the following system

[TABLE]

imply an optimal solution $\tilde{\boldsymbol{\upnu}}$ with $\tilde{\upnu}_{0}=0$ if and only if the vector of OLS residuals $\boldsymbol{e}$ belongs to the column space of $\mathrm{V}$. Since probability of $\boldsymbol{e}\in\mathscr{L}(\mathrm{V})$ is zero, our existence conclusions can be summarized in the next proposition.

Proposition 4 (Existence of NN-(M)DOOLSE)

In an orthogonal FDSLRM the following holds

1. $(0)$

NN-(M)DOOLSE* problems are strictly quadratic optimization problems.* 2. $(0)$

Their solutions $\tilde{\boldsymbol{\upnu}}$ always exist and they are unique global minimizers. 3. $(0)$

$\tilde{\upnu}_{0}=0\Leftrightarrow\boldsymbol{e}\in\mathscr{L}(\mathrm{V})\Leftrightarrow\boldsymbol{x}\in\mathscr{L}(\mathrm{F},\mathrm{V})$. 4. $(0)$

$P(\tilde{\upnu}_{0}=0)=P(\boldsymbol{e}\in\mathscr{L}(\mathrm{V}))=P\left(\boldsymbol{x}\in\mathscr{L}(\mathrm{F},\mathrm{V})\right)=0$.

KKT algorithm for NN-(M)DOOLSE

If we consider $\upnu_{0}=0$ which occurs if and only if $\boldsymbol{e}=\sum\limits_{j=1}^{l}\alpha_{j}\boldsymbol{v}_{j}$ in KKT conditions (6.38)-(6.40) then the optimal solution $\tilde{\boldsymbol{\upnu}}_{+}$ for NN-(M)DOOLSE is trivial: $\upnu_{+0}=0,\upnu_{+j}=\alpha_{j}^{2},j=1,\ldots,l$.

For $\upnu_{0}\neq 0$ implying $\uplambda_{0}=0$, $\upnu_{j}$ and $\uplambda_{j}$ cannot be simultaneous zero otherwise it would lead to contradiction between (6.38) and (6.39). Therefore, in the case $\upnu_{0}\neq 0$, the complementary slackness condition (3) $\upnu_{j}\uplambda_{j}=0,\,\,j\in\{1,\ldots,l\}$ can be rewritten in the form $\,b_{j}\upnu_{j}(1-b_{j})\uplambda_{j}=0$, where $b_{j}$ is an auxiliary indicator, which is zero if $\upnu_{j}=0,\uplambda_{j}\neq 0$ and one if $\upnu_{j}\neq 0,\uplambda_{j}=0$.

Using vector $\boldsymbol{b}=(b_{1},b_{2},\ldots,b_{l})^{\prime}\in\{0,1\}^{l}$, the derived KKT conditions (2) in Proposition 3 can be described by a $(l+1)\times(l+1)$ matrix function $\mathrm{K}(\boldsymbol{b})=\{K_{ij}(\boldsymbol{b})\}$ and a vector function $\boldsymbol{\upgamma}(\boldsymbol{b})=(\upgamma_{0}(\boldsymbol{b}),\upgamma_{1}(\boldsymbol{b}),\ldots,\upgamma_{l}(\boldsymbol{b}))^{\prime}$ as

[TABLE]

where

[TABLE]

Applying the Banachiewicz formula, we can write for the inverse of $\mathrm{K}(\boldsymbol{b})$ the following analytic expression

[TABLE]

where

$\mathrm{D}_{\boldsymbol{b}}=\mathrm{diag}\left\{b_{j}\left\|\boldsymbol{v}_{j}\right\|^{4}+b_{j}-1\right\}$, $\mathrm{G}_{\mathrm{V}}=\mathrm{diag}\left\{\left\|\boldsymbol{v}_{j}\right\|^{2}\right\}$,

$\boldsymbol{j}_{l}=(1,1,\ldots,1)^{\prime}\in\mathds{R}^{l}$, $\phi=n^{*}-\boldsymbol{b}^{\prime}\mathrm{G}_{\mathrm{V}}\mathrm{D}_{\boldsymbol{b}}^{-1}\mathrm{G}_{\mathrm{V}}\boldsymbol{j}_{l}$.

Finally, Proposition 4 guarantees the existence of unique auxiliary vectors ${\tilde{\boldsymbol{\upgamma}}},{\tilde{\boldsymbol{b}}}$:

[TABLE]

[TABLE]

Based on vectors ${\tilde{\boldsymbol{\upgamma}}},{\tilde{\boldsymbol{b}}}$, the NN-(M)DOOLSE $\tilde{\boldsymbol{\upnu}}_{+}=({\tilde{\upnu}}_{+0},{\tilde{\upnu}}_{+1},\ldots,{\tilde{\upnu}}_{+l})^{\prime}$ of $\boldsymbol{\upnu}$ as a solution of KKT conditions has the final form

[TABLE]

Thanks to Proposition 4, it is worth to mention that the matrix system (6.41) includes also solutions with $\upnu_{0}=0$.

Maximum likelihood estimators (RE)MLE

Generally, (RE)MLE in FDSLRM are not convex optimization problems with respect to $\boldsymbol{\upnu}$ and they do not exist for $\upnu_{0}=0$ which occurs if and only if $\boldsymbol{e}\in\mathscr{L}(V)$. If we apply the following bijective transformation defined on $\Upupsilon=(0,\infty)\times[0,\infty)^{l}$

[TABLE]

we can convert the (RE)MLE problems in the form of equivalent convex problems whose solutions can be readily converted to a (RE)MLE solutions by the inverse transformation to (6.42). Here are more detailed steps of the conversion.

In orthogonal FDSLRM $(\mathrm{F}^{\prime}\mathrm{V}=0,\mathrm{V}^{\prime}\mathrm{V}=\mathrm{diag}\left\{\left\|\boldsymbol{v}_{j}\right\|^{2}\right\})$, $\mathrm{\Sigma_{\boldsymbol{\upnu}}^{-1}}$ (stulajter_mse_2003, , lem 2.1) and $\boldsymbol{e}^{\prime}\mathrm{\Sigma_{\boldsymbol{\upnu}}^{-1}}\boldsymbol{e}$ are equal to

[TABLE]

Using orthogonality conditions (3.9) and expression (6.33), we get for determinants in (RE)ML loglikelihoods (4.25), (4.29)

[TABLE]

Since the chosen bijective transformation (6.42) also transforms convex constraints for $\boldsymbol{\upnu}$ to convex constraints for $\boldsymbol{d}$, substituting (6.43), (6.44) into objective functions (4.25), (4.29) of (RE)MLE, we can formulate the following proposition.

Proposition 5 (Equivalent (RE)MLE convex problems)

Let assume $\boldsymbol{e}\not\in\mathscr{L}(\mathrm{V})$ and consider a bijective transformation in the following form:

[TABLE]

Then, in a Gaussian orthogonal FDSLRM the (RE)MLE problems are equivalent to convex problems:

[TABLE]

where for $\emph{MLE: }n^{*}\equiv n,$ and $\emph{for REMLE: }n^{*}\equiv n-k$.

Once again, we can write the equivalent (RE)MLE problem with $2l$ constraints in the corresponding Lagrangian multiplier form based on the following Lagrangian with $2l$ multipliers $(\lambda_{1},\ldots,\lambda_{l},\mu_{1},\ldots,\mu_{l})^{\prime}=(\boldsymbol{\uplambda}^{\prime},\boldsymbol{\mu}^{\prime})^{\prime}\in\mathds{R}^{2l}_{+}$

[TABLE]

Then by the direct computation, we obtain KKT conditions for equivalent (RE)MLE describing (1) primal and dual feasibility for $\boldsymbol{d}$ and $(\boldsymbol{\uplambda}^{\prime},\boldsymbol{\mu}^{\prime})^{\prime}$, (2) stationarity of the Lagrangian ($\nabla_{\boldsymbol{d}}L(\boldsymbol{d},\boldsymbol{\uplambda},\boldsymbol{\mu})=0$) and (3) complementary slackness, all described by the next proposition.

Proposition 6 (KKT optimality conditions for equivalent (RE)MLE)

Consider equivalent convex optimization problems to (RE)MLE in the Lagrangian multiplier form

[TABLE]

Then, a necessary and sufficient condition for $\hat{\boldsymbol{d}},\hat{\boldsymbol{\uplambda}},\hat{\boldsymbol{\mu}}$ to be problems’ optimal solution is

1. $(0)$

$d_{0}-d_{j}\left\|\boldsymbol{v}_{j}\right\|^{2}>0,\,\,d_{j}\geq 0,\,\,\lambda_{j}\geq 0,\,\,\mu_{j}\geq 0,\,\,j\in\{1,\ldots,l\}$** 2. $(0)$

$||\boldsymbol{e}||^{2}-\dfrac{n^{*}-l}{d_{0}}-\displaystyle\sum\limits_{j=1}^{l}\left(\mu_{j}+\dfrac{1}{d_{0}-d_{j}||\boldsymbol{v}_{j}||^{2}}\right)=0,$* *

$\qquad\dfrac{||\boldsymbol{v}_{j}||^{2}}{d_{0}-d_{j}||\boldsymbol{v}_{j}||^{2}}-(\boldsymbol{e}^{\prime}\boldsymbol{v}_{j})^{2}-\lambda_{j}+\mu_{j}||\boldsymbol{v}_{j}||^{2}=0,\,\,j\in\{1,\ldots,l\}$** 3. $(0)$

$-d_{j}\lambda_{j}=0,\,\,-(d_{0}-d_{j}||\boldsymbol{v}_{j}||^{2})\mu_{j}=0,j\in\{1,\ldots,l\}$.

Similarly as in the case of NN-(M)-DOOLSE, the Hessian $\mathrm{H}=\nabla_{\boldsymbol{d}}^{2}f_{0}(\boldsymbol{d})$ equal to

[TABLE]

leads to strict convexity of the problem. Since $f_{0}(\boldsymbol{d})$ is also coercive, we can summarize the existence conditions of optimal solutions in the final proposition as we did for NN-(M)DOOLSE problems.

Proposition 7 (Existence of equivalent (RE)MLE)

Let assume $\boldsymbol{e}\not\in\mathscr{L}(\mathrm{V})$. Then in a Gaussian orthogonal FDSLRM the following holds

1. $(0)$

Equivalent (RE)MLE problems are strictly convex optimization problems. 2. $(0)$

The objective function $f_{0}(\boldsymbol{d})$ is coercive with respect to constraints. 3. $(0)$

Their solutions $\hat{\boldsymbol{d}}$ always exist and they are unique global minimizers.

Proof (Theorem 4.1)

The bijection (6.42) implies $d_{0}-d_{j}\left\|\boldsymbol{v}_{j}\right\|^{2}=(\upnu_{0}+\upnu_{j}\left\|\boldsymbol{v}_{j}\right\|^{2})^{-1}\neq 0$ which dictates $\mu_{j}=0$ to satisfy the complementary slackness conditions in KKT (3) and simultaneously forms the new version of KKT (1-3) of Proposition 6

[TABLE]

But for $\boldsymbol{e}\not\in\mathscr{L}(\mathrm{V})$, which occurs with probability 1, the system (6.49)-(6.51) is the same as the system (6.38)-(6.40), which is what we set out to prove.