  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.
