Authors: Jozef Hanč, Martina Hančová, Andrej Gajdoš
Faculty of Science, P. J. Šafárik University in Košice, Slovakia
emails: martina.hancova@upjs.sk
EBLUP-NE for electricity consumption 2
Python-based computational tools - SageMath
To get back to the contents, use the Home key.
SageMath is a free Python-based mathematics software, an open source alternative to the well-known commercial computer algebra systems Mathematica or Maple, which is built on top of SciPy ecosystem, Maxima, R, GAP, FLINT and many others open-source packages; http://www.sagemath.org/
# arbitrary precision in dps decimal places for numerical results
dps = 20
# precision in bits, real (floating-point) numbers with given precision
bits = round((dps+1)*ln(10)/ln(2))
RRk = RealField(bits)
# simplifying any symbolic matrix, vector
def simplify_matrix(A): return A.apply_map(lambda item: item.expand())
# converting symbolic matrix, vector to numerical in RRk
def RRk_matrix(A): return A.apply_map(RRk)
from itertools import product
In this FDSLRM application, we model the econometric time series data set, representing the hours observations of the consumption of the electric energy in some department store. The number of time series observations is $n=24$. The data and model were adapted from Štulajter & Witkovský, 2004.
The consumption data can be fitted by the following Gaussian orthogonal FDSLRM:
$ \begin{array}{rl} & X(t) & \! = \! & \beta_1+\beta_2\cos\left(\tfrac{2\pi t}{24}\right)+\beta_3\sin\left(\tfrac{2\pi t}{24}\right) +\\ & & & +Y_1\cos\left(\tfrac{2\pi t\cdot 2}{24}\right)+Y_2\sin\left(\tfrac{2\pi t\cdot 2}{24}\right)+\\ & & & +Y_3\cos\left(\tfrac{2\pi t\cdot 3}{24}\right)+Y_4\sin\left(\tfrac{2\pi t\cdot 3}{24}\right) +w(t), \, t\in \mathbb{N}, \end{array} $
where $\boldsymbol{\beta}=(\beta_1,\,\beta_2,\,\beta_3)' \in \mathbb{R}^3\,,\mathbf{Y} = (Y_1, Y_2, Y_3, Y_4)' \sim \mathcal{N}_4(\boldsymbol{0}, \mathrm{D})\,, w(t) \sim \mathcal{iid}\, \mathcal{N} (0, \sigma_0^2)\,, \boldsymbol{\nu}= (\sigma_0^2, \sigma_1^2, \sigma_2^2, \sigma_3^2, \sigma_4^2) \in \mathbb{R}_{+}^5.$
# data - time series observation x
data = [40.3,40.7,38.5,37.9,38.6,41.1,45.2,45.7,46.7,46.5,
45.2,45.1,45.8,46.3,47.5,48.5,49.1,51.7,50.6,48.0,
44.7,41.2,40.0,40.3]
# observation x as a one-column matrix of rational numbers (infinite precision)
x = matrix(QQ, data).T
xc = x.T
xc
[403/10 407/10 77/2 379/10 193/5 411/10 226/5 457/10 467/10 93/2 226/5 451/10 229/5 463/10 95/2 97/2 491/10 517/10 253/5 48 447/10 206/5 40 403/10]
# model parameters
n, k, l = 24, 3, 4
# significant frequencies
om1, om2, om3 = 2*pi/24, 2*pi/12, 2*pi/8
# model design matrices F', F, V',V
Fc = matrix([[1 for t in range(1,n+1)],
[cos(om1*t) for t in range(1,n+1)],
[sin(om1*t) for t in range(1,n+1)]])
Vc = matrix([[cos(om2*t) for t in range(1,n+1)],
[sin(om2*t) for t in range(1,n+1)],
[cos(om3*t) for t in range(1,n+1)],
[sin(om3*t) for t in range(1,n+1)]])
F, V = Fc.T, Vc.T
# columns vj of V and their squared norm ||vj||^2
vv = lambda j: V[:,j-1]
nv2 = lambda j: (vv(j).T*vv(j)).trace().expand()
# auxiliary matrices and vectors
# Gram matrices GF, GV
GF, GV = (Fc*F).expand(), (Vc*V).expand()
InvGF, InvGV = (GF^(-1)).expand(), (GV^(-1)).expand()
# projectors PF, MF, PV, MV
In = identity_matrix(n)
PF = (F*InvGF*Fc).canonicalize_radical().expand()
PV = (V*InvGV*Vc).canonicalize_radical().expand()
MF, MV = (In-PF).expand(), (In-PV).expand()
# residuals e, e'
e = MF*x
ec = e.T
# orthogonality condition
show((Fc*V).expand())
show(GV)
using formula (4.1) from Hancova et al 2019
$
\renewcommand{\arraystretch}{1.4} \breve{\boldsymbol{\nu}}(\mathbf{e}) = \begin{pmatrix} \tfrac{1}{n-k-l}\,\mathbf{e}'\,\mathrm{M_V}\,\mathbf{e} \\ (\mathbf{e}'\mathbf{v}_1)^2/||\mathbf{v}_1||^4 \\ \vdots \\ (\mathbf{e}'\mathbf{v}_l)^2/||\mathbf{v}_l||^4 \end{pmatrix} $
# NE according to formula (4.1)
NE0 = [1/(n-k-l)*(ec*MV*e).trace()]
NEj = [((ec*vv(j)).trace()/nv2(j))^2 for j in [1..l]]
NE = vector(NE0+NEj)
# final NE
NEsimp = NE.column().expand()
show(NEsimp)
# numerical results with given precision in dps
NEnum = RRk_matrix(NEsimp)
NEnum
[ 1.0930446920400417197] [ 2.9657173646433129174] [ 1.7618587371177719801] [0.37193497450591316960] [ 1.8634794260764497182]
NEnum.norm()
4.087207309639791
using formula (3.10) from Hancova et al 2019.
$
\mathring{\nu}_j = \rho_j^2 \breve{\nu}_j; j = 0,1 \ldots, l\ \rho_0 = 1, \rho_j = \dfrac{\hat{\nu}_j||\mathbf{v}_j||^2}{\hat{\nu}_0+\hat{\nu}_j||\mathbf{v}_j||^2} $
where $\boldsymbol{\breve{\nu}}$ are NE, $\boldsymbol{\hat{\nu}}$ are initial estimates for EBLUP-NE
# EBLUP-NE based on formula (3.9)
rho2 = lambda est: vector( [1] + [ (est[j]*nv2(j)/(est[0]+est[j]*nv2(j)))^2 for j in [1..l] ])
EBLUPNE = lambda est: vector( rho2(est).pairwise_product(NE).list() )
# final EBLUPNE based on NE
EBLUPNEsimp = EBLUPNE(NE).column().expand()
show(EBLUPNEsimp)
# numerical results
rho2(NEnum).column()
[ 1] [0.94129166721252861086] [0.90410056203431227372] [0.64525402743480536407] [0.90896740457337881883]
EBLUPNEnum = RRk_matrix(EBLUPNEsimp)
EBLUPNEnum
[ 1.0930446920400417197] [ 2.7916050426462506682] [ 1.5928974744532412866] [0.23999254024380213000] [ 1.6938420573966000382]
EBLUPNEnum.norm()
3.801555617352273
using the the KKT algorithm (tab.3, Hancova et al 2019)
$~$
$
\qquad \mathbf{q} = \left(\begin{array}{c} \mathbf{e}' \mathbf{e}\\ (\mathbf{e}' \mathbf{v}_{1})^2 \\ \vdots \\ (\mathbf{e}' \mathbf{v}_{l})^2 \end{array}\right) $
$\qquad\mathrm{G} = \left(\begin{array}{ccccc}
\small n^* & ||\mathbf{v}{1}||^2 & ||\mathbf{v}{2}||^2 & \ldots & ||\mathbf{v}{l}||^2 \ ||\mathbf{v}{1}||^2 & ||\mathbf{v}{1}||^4 & 0 & \ldots & 0 \ ||\mathbf{v}{2}||^2 & 0 & ||\mathbf{v}{2}||^4 & \ldots & 0 \ \vdots & \vdots & \vdots & \ldots & \vdots \ ||\mathbf{v}{l}||^2 & 0 & 0 & \ldots & ||\mathbf{v}_{l}||^4 \end{array}\right) $
# Input: form G
ns = n
u, v, Q = ns, matrix([nv2(j) for j in [1..l]]), diagonal_matrix([nv2(j)^2 for j in [1..l]])
G = block_matrix([[u,v],[v.T,Q]])
# form q
e2, Ve2 = ec*e, (Vc*e).elementwise_product(Vc*e)
q = e2.stack(Ve2).expand()
# body of the algorithm
for b in product([0,1], repeat=l+1):
# set the KKT-conditions matrix K
K = G*1
for j in range(l+1):
if b[j] == 1: K[0,j], K[j,j] = 0,-1
# calculate the auxiliary vector g
g = K^(-1)*q
# test non-negativity g
if RRk(min(g.list()))>=0: break
RRk_matrix(K^(-1)*q), b
( [0.92908798823403546177] [ 2.8882933656238099623] [ 1.6844347380982690249] [0.29451097548641021445] [ 1.7860554270569467631], (0, 0, 0, 0, 0) )
# Output: Form estimates nu
nu = g*1
for j in range(l+1):
if b[j] == 1: nu[j] = 0
# final DOOLSE
DOOLSE = nu
show(DOOLSE)
# numerical results
DOOLSEnum = RRk_matrix(DOOLSE)
print(DOOLSEnum, DOOLSE.norm(),b)
([0.92908798823403546177] [ 2.8882933656238099623] [ 1.6844347380982690249] [0.29451097548641021445] [ 1.7860554270569467631], 3.9140125377802497, (0, 0, 0, 0, 0))
# final EBLUPNE based on DOOLSE
EBLUPNEsimp = EBLUPNE(DOOLSE).column().expand()
show(EBLUPNEsimp)
# numerical results
rho2(DOOLSEnum).column()
[ 1] [0.94846887861963596617] [0.91404212152622161047] [0.62700200702811399302] [0.91863007601174459834]
EBLUPNEnum = RRk_matrix(EBLUPNEsimp)
EBLUPNEnum
[ 1.0930446920400417197] [ 2.8128906231460250176] [ 1.6104130979046378695] [0.23320397549915796799] [ 1.7118482468229312038]
EBLUPNEnum.norm()
3.832145510914468
using the KKT algorithm (tab.3, Hancova et al 2019)
# Input: form G
ns = n-k
u, v, Q = ns, matrix([nv2(j) for j in [1..l]]), diagonal_matrix([nv2(j)^2 for j in [1..l]])
G = block_matrix([[u,v],[v.T,Q]])
# form q
e2, Ve2 = ec*e, (Vc*e).elementwise_product(Vc*e)
q = e2.stack(Ve2).expand()
# body of the algorithm
for b in product([0,1], repeat=l+1):
# set the KKT-conditions matrix K
K = G*1
for j in range(l+1):
if b[j] == 1: K[0,j], K[j,j] = 0,-1
# calculate the auxiliary vector g
g = K^(-1)*q
# test non-negativity g
if RRk(min(g.list()))>=0: break
# Output: Form estimates nu
nu = g*1
for j in range(l+1):
if b[j] == 1: nu[j] = 0
# final MDOOLSE
MDOOLSE = nu
show(MDOOLSE)
# numerical results
MDOOLSEnum = RRk_matrix(nu)
print(MDOOLSEnum, MDOOLSE.norm(),b)
([ 1.0930446920400417197] [ 2.8746303069733094408] [ 1.6707716794477685035] [0.28084791683590969296] [ 1.7723923684064462416], 3.933188829103883, (0, 0, 0, 0, 0))
# final EBLUPNE based on NE
EBLUPNEsimp = EBLUPNE(MDOOLSE).column().expand()
show(EBLUPNEsimp)
# numerical results
rho2(MDOOLSEnum).column()
[ 1] [0.93951664761802805166] [0.89927400850167713940] [0.57017526933078310576] [0.90462906726347676354]
EBLUPNEnum = RRk_matrix(EBLUPNEsimp)
EBLUPNEnum
[ 1.0930446920400417197] [ 2.7863408362122582278] [ 1.5843937689416014278] [0.21206812426244698957] [ 1.6857576550762177074]
EBLUPNEnum.norm()
3.78886491318682
This notebook belongs to suplementary materials of the paper submitted to Statistical Papers and available at https://arxiv.org/abs/1905.07771.
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.
Discrete Spectrum Linear Regression Models](https://link.springer.com/article/10.1007/s001840300299). Metrika, 2004, Vol. 60, No. 2, pp. 105–118