PY-estimation-electricity2-SciPyCVXPY.ipynb

Kernel: Python 2 (Ubuntu Linux)

| Table of Contents | Data and model | Natural estimators | NN-DOOLSE, MLE | NN-MDOOLSE, REMLE | References |

Authors: Jozef Hanč, Martina Hančová, Andrej Gajdoš
Faculty of Science, P. J. Šafárik University in Košice, Slovakia
emails: [email protected]

---

EBLUP-NE for electricity consumption 2

Python-based computational tools - SciPy, CVXPY

Table of Contents

• Data and model - data and model description of empirical data

• Natural estimators - EBLUPNE based on NE

• NN-DOOLSE, MLE - EBLUPNE based on nonnegative DOOLSE (same as MLE)

• NN-MDOOLSE, REMLE - EBLUPNE based on nonnegative MDOOLSE (same as REMLE)

• References

To get back to the contents, use the Home key.

Python libraries

CVXPY: A Python-Embedded Modeling Language for Convex Optimization

• Purpose: scientific Python library for solving convex optimization tasks

• Version: 1.0.1, 2018

• URL: https://www.cvxpy.org/

• Computational parameters of CVXPY:

• solver - the convex optimization solver ECOS, OSQP, and SCS chosen according to the given problem * OSQP for convex quadratic problems * max_iter - maximum number of iterations (default: 10000). * eps_abs - absolute accuracy (default: 1e-4). * eps_rel - relative accuracy (default: 1e-4). * ECOS for convex second-order cone problems * max_iters - maximum number of iterations (default: 100). * abstol - absolute accuracy (default: 1e-7). * reltol - relative accuracy (default: 1e-6). * feastol - tolerance for feasibility conditions (default: 1e-7). * abstol_inacc - absolute accuracy for inaccurate solution (default: 5e-5). * reltol_inacc - relative accuracy for inaccurate solution (default: 5e-5). * feastol_inacc - tolerance for feasibility condition for inaccurate solution (default: 1e-4). * SCS for large-scale convex cone problems * max_iters - maximum number of iterations (default: 2500). * eps - convergence tolerance (default: 1e-4). * alpha - relaxation parameter (default: 1.8). * scale - balance between minimizing primal and dual residual (default: 5.0). * normalize - whether to precondition data matrices (default: True). * use_indirect - whether to use indirect solver for KKT sytem (instead of direct) (default: True).

Scipy - NumPy, Pandas

• Numpy is the fundamental Python library of SciPy ecosystem for fast scientific computing with large, multi-dimensional arrays and matrices, along with a large collection of high-level mathematical functions to operate on these arrays.

• default precision: double floating-point precision $\text{eps}<10^{-15}$

• Pandas is the Python library providing high-performance, easy to use data structures.

import cvxpy
import numpy as np
import pandas as pd
import platform as pt

from cvxpy import *
from math import cos, sin
from numpy.linalg import inv, norm
from itertools import product
from __future__ import print_function
from __future__ import division

np.set_printoptions(precision=10)

# software versions
print('cvxpy:'+cvxpy.__version__)
print('numpy:'+np.__version__)
print('pandas:'+pd.__version__)
print('python:'+pt.python_version())

cvxpy:1.0.22
numpy:1.16.3
pandas:0.23.4
python:2.7.15rc1

---

Data and Toy Model 2

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.$

SciPy(Numpy)

# data - time series observation
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 matrix
x = np.mat(data).T
xc = x.T

# model parameters
n, k, l = 24, 3, 4

# significant frequencies
om1, om2, om3 = 2*np.pi/24, 2*np.pi/12, 2*np.pi/8

# model - design matrices F', F, V',V
Fc = np.mat([[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 = np.mat([[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: np.trace(vv(j).T*vv(j))

# auxiliary matrices and vectors

# Gram matrices GF, GV
GF, GV = Fc*F, Vc*V
InvGF, InvGV = inv(GF), inv(GV)

# projectors PF, MF, PV, MV
In = np.identity(n)
PF = F*InvGF*Fc
PV = V*InvGV*Vc
MF, MV = In-PF, In-PV

# residuals e, e'
e = MF*x
ec = e.T

# orthogonality condition
Fc*V, GV

(matrix([[-1.7763568394e-15,  3.9831970004e-16, -2.2204460493e-15,
          -1.2898995918e-15],
         [-2.4424906542e-15, -1.3843693787e-15, -3.1086244690e-15,
          -5.1870949789e-16],
         [-7.2154905439e-16, -2.2782998587e-16,  1.0883751302e-15,
           4.3896384668e-16]]),
 matrix([[ 1.2000000000e+01,  4.8066158344e-16,  2.2204460493e-16,
          -2.9836905420e-15],
         [ 4.8066158344e-16,  1.2000000000e+01, -8.4428500303e-16,
           1.2981158543e-15],
         [ 2.2204460493e-16, -8.4428500303e-16,  1.2000000000e+01,
           1.7178491110e-15],
         [-2.9836905420e-15,  1.2981158543e-15,  1.7178491110e-15,
           1.2000000000e+01]]))

---

Natural estimators

ANALYTICALLY

using formula (4.1) from Hancova et al 2019

ParseError: KaTeX parse error: \renewcommand{\arraystretch} when command \arraystretch does not yet exist; use \newcommand

$\boldsymbol{1^{st}}$ stage of EBLUP-NE

SciPy(Numpy)

# NE according to formula (4.1)
NE0 = [1/(n-k-l)*np.trace(ec*MV*e)]
NEj = [(np.trace(ec*vv(j))/nv2(j))**2 for j in range(1,l+1)] 
NE = NE0+NEj
print(NE, norm(NE))

[1.093044692040042, 2.965717364643324, 1.7618587371177696, 0.3719349745059126, 1.8634794260764524] 4.0872073096398

CVXPY

NE as a convex optimization problem

$\begin{array}{ll} \textit{minimize} & \quad f_0(\boldsymbol{\nu})=||\mathbf{e}\mathbf{e}' - \mathrm{VDV'}||^2+||\mathrm{M_V}\mathbf{e}\mathbf{e}'\mathrm{M_V}-\nu_0\mathrm{M_F}\mathrm{M_V}||^2 $6pt] \textit{subject to} & \quad \boldsymbol{\nu} = \left(\nu_0, \ldots, \nu_l\right)'\in [0, \infty)^{l+1} \end{array}$

# the optimization variable, objective function
v = Variable(l+1)
fv = sum_squares(e*ec-V*diag(v[1:])*Vc)+sum_squares(MV*e*ec*MV-v[0]*MF*MV)

# the optimization problem for NE
objective = Minimize(fv)
constraints = [v >= 0]
prob = Problem(objective,constraints)

# solve the NE problem
prob.solve()
NEcvx = v.value
print('NEcvx =', NEcvx, ' norm = ', norm(NEcvx))

NEcvx = [1.093044692  2.9657173646 1.7618587371 0.3719349745 1.8634794261]  norm =  4.087207309639801

$\boldsymbol{2^{nd}}$ stage of EBLUP-NE

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

SciPy(Numpy)

# EBLUP-NE based on formula (3.9)
rho2 = lambda est: [1]+ [ (est[j]*nv2(j)/(est[0]+est[j]*nv2(j)))**2 for j in range(1,l+1) ]
EBLUPNE = lambda est: [rho2(est)[j]*NE[j] for j in range(l+1)]

# numerical results
print(rho2(NE))

[1, 0.9412916672125291, 0.9041005620343122, 0.6452540274348048, 0.908967404573379]

print(EBLUPNE(NE),norm(EBLUPNE(NE)))

[1.093044692040042, 2.7916050426462626, 1.5928974744532391, 0.23999254024380154, 1.6938420573966029] 3.8015556173522818

cross-checking

using formula (3.6) for general FDSLRM from Hancova et al 2019.

$\mathring{\nu}_0 = \breve{\nu}_0, \mathring{\nu}_j = (\mathbf{Y}^*)_j^2, j = 1, 2, \ldots, l \\ \mathbf{Y}^* = \mathbb{T}\mathbf{X} \mbox{ with } \mathbb{T} = \mathrm{D}\mathbb{U}^{-1}\mathrm{V}'\mathrm{M_F}, \mathbb{U} = \mathrm{V}'\mathrm{M_F}\mathrm{V}\mathrm{D} + \nu_0 \mathrm{I}_l$

def EBLUPNEgen(est):
    D = np.diag(est[1:])
    U = Vc*MF*V*D + est[0]*np.identity(l)
    T = D*inv(U)*Vc*MF
    eest = np.vstack((np.matrix(NE[0]),np.multiply(T*x, T*x)))
    return np.array(eest).flatten().tolist()

print(EBLUPNEgen(NE), norm(EBLUPNEgen(NE)))

[1.093044692040042, 2.7916050426462733, 1.5928974744532476, 0.23999254024379996, 1.6938420573965995] 3.8015556173522915

---

NN-DOOLSE or MLE

$\boldsymbol{1^{st}}$ stage of EBLUP-NE

KKT algorithm

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)$

SciPy(Numpy)

## SciPy(Numpy)# Input: form G
ns, nvj = n, norm(V, axis=0)
u, v, Q  = np.mat(ns), np.mat(nvj**2), np.diag(nvj**4)
G = np.bmat([[u,v],[v.T,Q]])
# form q
e2, Ve2 = ec*e, np.multiply(Vc*e, Vc*e)
q = np.vstack((e2, Ve2))

# body of the algorithm
for b in product([0,1], repeat=l): 
    # set the KKT-conditions matrix K
    K = G*1
    for j in range(1,l+1): 
        if b[j-1] == 0: K[0,j], K[j,j]  = 0,-1
    # calculate the auxiliary vector g
    g = inv(K)*q
    # test non-negativity g
    if (g >= 0).all(): break   

# Output: Form estimates nu
nu = g*1
for j in range(1,l+1):
    if b[j-1] == 0: nu[j] = 0

NN_DOOLSE = np.array(nu).flatten()
print(NN_DOOLSE, norm(NN_DOOLSE),b)

[0.9290879882 2.8882933656 1.6844347381 0.2945109755 1.7860554271] 3.91401253778026 (1, 1, 1, 1)

CVXPY

nonnegative DOOLSE as a convex optimization problem

ParseError: KaTeX parse error: Got function '\boldsymbol' with no arguments as subscript at position 97: …hbf{e}'-\Sigma_\̲b̲o̲l̲d̲s̲y̲m̲b̲o̲l̲{\nu}||^2 $6pt]…

# set the optimization variable, objective function
v = Variable(l+1)
fv = sum_squares(e*e.T - v[0]*In - V*diag(v[1:])*Vc)

# construct the problem for DOOLSE
objective = Minimize(fv)
constraints = [v >= 0]
prob = Problem(objective,constraints)

# solve the DOOLSE problem
prob.solve()
NN_DOOLSEcvx = v.value
print('NN-DOOLSEcvx =', NN_DOOLSEcvx, 'norm =', norm(NN_DOOLSEcvx))

NN-DOOLSEcvx = [0.9290879882 2.8882933656 1.6844347381 0.2945109755 1.7860554271] norm = 3.9140125377802586

CVXPY

using equivalent (RE)MLE convex problem (proposition 5, Hancova et al 2019)

$\begin{array}{ll} \textit{minimize} & \quad f_0(\mathbf{d})=-(n^*\!-l)\ln d_0 - \displaystyle\sum\limits_{j=1}^{l} \ln(d_0-d_j||\mathbf{v}_j||^2+d_0\mathbf{e}'\mathbf{e}-\mathbf{e}'\mathrm{V}\,\mathrm{diag}\{d_j\}\mathrm{V}'\mathbf{e} $6pt] \textit{subject to} & \quad d_0 > \max\{d_j||\mathbf{v}_j||^2, j = 1, \ldots, l\} \\ & \quad d_j \geq 0, j=1,\ldots l \\ & \\ & \quad\text{for MLE: } n^* = n, \text{ for REMLE: } n^* = n-k \\ \textit{back transformation:} & \quad \nu_0 = \dfrac{1}{d_0}, \nu_j = \dfrac{d_j}{d_0\left(d_0 -d_j||\mathbf{v}_j||^2\right)} \end{array}$

# set variables for the objective
ns = n
d = Variable(l+1)
logdetS = (ns-l)*log(d[0])+sum(log(d[0]-GV*d[1:]))

# construct the problem
objective = Maximize(logdetS-(d[0]*ec*e-ec*V*diag(d[1:])*Vc*e))
constraints = [0 <= d[1:], max(GV*d[1:]) <= d[0]]
prob = Problem(objective,constraints)

# solve the problem
solution = prob.solve()
dv = d.value.tolist()

# back transformation
s0 = [1/dv[0]]
sj = [dv[i]/(dv[0]*(dv[0]-dv[i]*GV[i-1,i-1])) for i in range(1,l+1)]
sv = s0+sj

print('REMLEcvx = ', sv, ' norm = ', norm(sv))

REMLEcvx =  [0.9290880038928567, 2.8882936875435075, 1.6844348208009916, 0.2945109888857377, 1.7860555721680658]  norm =  3.914012881871393

$\boldsymbol{2^{nd}}$ stage of EBLUP-NE

SciPy(Numpy)

## SciPy(Numpy)# numerical results
print(rho2(NN_DOOLSE))

[1, 0.948468878619636, 0.9140421215262213, 0.6270020070281136, 0.9186300760117446]

print(EBLUPNE(NN_DOOLSE), norm(EBLUPNE(NN_DOOLSE)))

[1.093044692040042, 2.812890623146036, 1.6104130979046352, 0.23320397549915747, 1.7118482468229337] 3.832145510914476

#cross-checking
print(EBLUPNEgen(NN_DOOLSE), norm(EBLUPNEgen(NN_DOOLSE)))

[1.093044692040042, 2.8128906231460324, 1.6104130979046318, 0.23320397549915772, 1.7118482468229368] 3.8321455109144735

---

NN-MDOOLSE or REMLE

using the KKT algorithm (tab.3, Hancova et al 2019)

$\boldsymbol{1^{st}}$ stage of EBLUP-NE

KKT algorithm

SciPy(Numpy)

# Input: form G
ns, nvj = n-k, norm(V, axis=0)
u, v, Q  = np.mat(ns), np.mat(nvj**2), np.diag(nvj**4)
G = np.bmat([[u,v],[v.T,Q]])
# form q
e2, Ve2 = ec*e, np.multiply(Vc*e, Vc*e)
q = np.vstack((e2, Ve2))

# body of the algorithm
for b in product([0,1], repeat=l): 
    # set the KKT-conditions matrix K
    K = G*1
    for j in range(1,l+1): 
        if b[j-1] == 0: K[0,j], K[j,j]  = 0,-1
    # calculate the auxiliary vector g
    g = inv(K)*q
    # test non-negativity g
    if (g >= 0).all(): break   

# Output: Form estimates nu
nu = g*1
for j in range(1,l+1):
    if b[j-1] == 0: nu[j] = 0

NN_MDOOLSE = np.array(nu).flatten()
print(NN_MDOOLSE, norm(NN_MDOOLSE),b)

[1.093044692  2.874630307  1.6707716794 0.2808479168 1.7723923684] 3.933188829103891 (1, 1, 1, 1)

CVXPY

nonnegative DOOLSE as a convex optimization problem

ParseError: KaTeX parse error: Got function '\boldsymbol' with no arguments as subscript at position 109: …hrm{M_F}\Sigma_\̲b̲o̲l̲d̲s̲y̲m̲b̲o̲l̲{\nu}\mathrm{M_…

# set the optimization variable, objective function
v = Variable(l+1)
fv = sum_squares(e*e.T - v[0]*MF - V*diag(v[1:])*V.T)

# construct the problem for DOOLSE
objective = Minimize(fv)
constraints = [v >= 0]
prob = Problem(objective,constraints)

# solve the DOOLSE problem
prob.solve()

NN_MDOOLSEcvx = v.value
print('NN_MDOOLSEcvx =', NN_MDOOLSEcvx, 'norm =', norm(NN_MDOOLSEcvx))

NN_MDOOLSEcvx = [1.093044692  2.874630307  1.6707716794 0.2808479168 1.7723923684] norm = 3.933188829103891

CVXPY

using equivalent (RE)MLE convex problem (proposition 5, Hancova et al 2019)

$\begin{array}{ll} \textit{minimize} & \quad f_0(\mathbf{d})=-(n^*\!-l)\ln d_0 - \displaystyle\sum\limits_{j=1}^{l} \ln(d_0-d_j||\mathbf{v}_j||^2+d_0\mathbf{e}'\mathbf{e}-\mathbf{e}'\mathrm{V}\,\mathrm{diag}\{d_j\}\mathrm{V}'\mathbf{e} $6pt] \textit{subject to} & \quad d_0 > \max\{d_j||\mathbf{v}_j||^2, j = 1, \ldots, l\} \\ & \quad d_j \geq 0, j=1,\ldots l \\ & \\ & \quad\text{for MLE: } n^* = n, \text{ for REMLE: } n^* = n-k \\ \textit{back transformation:} & \quad \nu_0 = \dfrac{1}{d_0}, \nu_j = \dfrac{d_j}{d_0\left(d_0 -d_j||\mathbf{v}_j||^2\right)} \end{array}$

# set variables for the objective
ns = n - k
d = Variable(l+1)
logdetS = (ns-l)*log(d[0])+sum(log(d[0]-GV*d[1:]))

# construct the problem
objective = Maximize(logdetS-(d[0]*ec*e-ec*V*diag(d[1:])*Vc*e))
constraints = [0 <= d[1:], max(GV*d[1:]) <= d[0]]
prob = Problem(objective,constraints)

# solve the problem
solution = prob.solve()
dv = d.value.tolist()

# back transformation
s0 = [1/dv[0]]
sj = [dv[i]/(dv[0]*(dv[0]-dv[i]*GV[i-1,i-1])) for i in range(1,l+1)]
sv = s0+sj

print('REMLEcvx = ', sv, ' norm = ', norm(sv))

REMLEcvx =  [1.0930447242738333, 2.874631059815253, 1.6707717649941174, 0.28084791305626483, 1.7723924736455892]  norm =  3.9331894717801954

$\boldsymbol{2^{nd}}$ stage of EBLUP-NE

SciPy(Numpy)

# numerical results
print(rho2(NN_MDOOLSE))

[1, 0.9395166476180281, 0.8992740085016769, 0.5701752693307821, 0.9046290672634769]

print(EBLUPNE(NN_MDOOLSE),norm(EBLUPNE(NN_MDOOLSE)))

[1.093044692040042, 2.786340836212269, 1.584393768941599, 0.2120681242624463, 1.6857576550762206] 3.788864913186828

#cross-checking
print(EBLUPNEgen(NN_MDOOLSE), norm(EBLUPNEgen(NN_MDOOLSE)))

[1.093044692040042, 2.7863408362122786, 1.5843937689416003, 0.21206812426244492, 1.685757655076218] 3.788864913186835

---

References

This notebook belongs to suplementary materials of the paper submitted to Statistical Papers and available at https://arxiv.org/abs/1905.07771.

• Hančová, M., Vozáriková, G., Gajdoš, A., Hanč, J. (2019). Estimating variance components in time series linear regression models using empirical BLUPs and convex optimization, https://arxiv.org/, 2019.

Abstract of the paper

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.

• Štulajter, F., Witkovský, V. (2004). Estimation of Variances in Orthogonal Finite Discrete Spectrum Linear Regression Models. Metrika, 2004, Vol. 60, No. 2, pp. 105–118

| Table of Contents | Data and model | Natural estimators | NN-DOOLSE, MLE | NN-MDOOLSE, REMLE | References |