Coverage for pybreakpoints/recresid_.py : 90%

# -*- coding: utf-8 -*- 

u""" Recursive residuals computation 

Citations: 

- Brown, RL, J Durbin, and JM Evans. 1975. Techniques for Testing the 

Consistency of Regression Relationships over Time. Journal of the Royal 

Statistical Society. Series B (Methodological) 37 (2): 149-192. 

- Judge George G., William E. Griffiths, R. Carter Hill, Helmut Lütkepohl, 

and Tsoung-Chao Lee. 1985. The theory and practice of econometrics. 

New York: Wiley. ISBN: 978-0-471-89530-5 

""" 

import numpy as np 

import pandas as pd 

import xarray as xr 

from .core import PANDAS_LIKE 

from .compat import jit 

@jit(nopython=True, nogil=True) 

def _recresid(X, y, span): 

nobs, nvars = X.shape 

recresid_ = np.nan * np.zeros((nobs)) 

recvar = np.nan * np.zeros((nobs)) 

X0 = X[:span, :] 

y0 = y[:span] 

# Initial fit 

XTX_j = np.linalg.inv(np.dot(X0.T, X0)) 

XTY = np.dot(X0.T, y0) 

beta = np.dot(XTX_j, XTY) 

yhat_j = np.dot(X[span - 1, :], beta) 

recresid_[span - 1] = y[span - 1] - yhat_j 

recvar[span - 1] = 1 + np.dot(X[span - 1, :], 

np.dot(XTX_j, X[span - 1, :])) 

for j in range(span, nobs): 

x_j = X[j:j+1, :] 

y_j = y[j] 

# Prediction with previous beta 

resid_j = y_j - np.dot(x_j, beta) 

# Update 

XTXx_j = np.dot(XTX_j, x_j.T) 

f_t = 1 + np.dot(x_j, XTXx_j) 

XTX_j = XTX_j - np.dot(XTXx_j, XTXx_j.T) / f_t # eqn 5.5.15 

beta = beta + (XTXx_j * resid_j / f_t).ravel() # eqn 5.5.14 

recresid_[j] = resid_j.item() 

recvar[j] = f_t.item() 

return recresid_ / np.sqrt(recvar) 

def recresid(X, y, span=None): 

""" Return standardized recursive residuals for y ~ X 

Parameters 

---------- 

X : array like 

2D (n_obs x n_features) design matrix 

y : array like 

1D independent variable 

span : int, optional 

Minimum number of observations for initial regression. If ``span`` 

is None, use the number of features in ``X`` 

Returns 

------- 

array like 

np.ndarray, pd.Series, or xr.DataArray containing recursive residuals 

standardized by prediction error variance 

Notes 

----- 

For a matrix :math:`X_t` of :math:`T` total observations of :math:`n` 

variables, the :math:`t` th recursive residual is the forecast prediction 

error for :math:`y_t` using a regression fit on the first :math:`t - 1` 

observations. Recursive residuals are scaled and standardized so they are 

:math:`N(0, 1)` distributed. 

Using notation from Brown, Durbin, and Evans (1975) and Judge, et al 

(1985): 

.. math:: 

w_r = 

\\frac{y_r - \\boldsymbol{x}_r^{\prime}\\boldsymbol{b}_{r-1}} 

{\sqrt{(1 + \\boldsymbol{x}_r^{\prime} 

S_{r-1}\\boldsymbol{x}_r)}} 

= 

\\frac 

{y_r - \\boldsymbol{x}_r^{\prime}\\boldsymbol{b}_r} 

{\sqrt{1 - \\boldsymbol{x}_r^{\prime}S_r\\boldsymbol{x}_r}} 

r = k + 1, \ldots, T, 

where :math:`S_{r}` is the residual sum of squares after 

fitting the model on :math:`r` observations. 

A quick way of calculating :math:`\\boldsymbol{b}_r` and 

:math:`S_r` is using an update formula (Equations 4 and 5 in 

Brown, Durbin, and Evans; Equation 5.5.14 and 5.5.15 in Judge et al): 

.. math:: 

\\boldsymbol{b}_r 

= 

b_{r-1} + 

\\frac 

{S_{r-1}\\boldsymbol{x}_j 

(y_r - \\boldsymbol{x}_r^{\prime}\\boldsymbol{b}_{r-1})} 

{1 + \\boldsymbol{x}_r^{\prime}S_{r-1}x_r} 

.. math:: 

S_r = 

S_{j-1} - 

\\frac{S_{j-1}\\boldsymbol{x}_r\\boldsymbol{x}_r^{\prime}S_{j-1}} 

{1 + \\boldsymbol{x}_r^{\prime}S_{j-1}\\boldsymbol{x}_r} 

See Also 

-------- 

statsmodels.stats.diagnostic.recursive_olsresiduals 

""" 

if not span: 

span = X.shape[1] 

_X = X.values if isinstance(X, PANDAS_LIKE) else X 

_y = y.values.ravel() if isinstance(y, PANDAS_LIKE) else y.ravel() 

rresid = _recresid(_X, _y, span)[span:] 

if isinstance(y, PANDAS_LIKE): 

if isinstance(y, (pd.Series, pd.DataFrame)): 

rresid = pd.Series(data=rresid, 

index=y.index[span:], 

name='recresid') 

elif isinstance(y, xr.DataArray): 

rresid = xr.DataArray(rresid, 

coords={'time': y.get_index('time')[span:]}, 

dims=('time', ), 

name='recresid') 

return rresid