code: def recursive_olsresiduals(
      res, skip=None, lamda=0.0, alpha=0.95, order_by=None
  ):
      """
      Calculate recursive ols with residuals and Cusum test statistic
  
      Parameters
      ----------
      res : RegressionResults
          Results from estimation of a regression model.
      skip : int, default None
          The number of observations to use for initial OLS, if None then skip is
          set equal to the number of regressors (columns in exog).
      lamda : float, default 0.0
          The weight for Ridge correction to initial (X'X)^{-1}.
      alpha : {0.90, 0.95, 0.99}, default 0.95
          Confidence level of test, currently only two values supported,
          used for confidence interval in cusum graph.
      order_by : array_like, default None
          Integer array specifying the order of the residuals. If not provided,
          the order of the residuals is not changed. If provided, must have
          the same number of observations as the endogenous variable.
  
      Returns
      -------
      rresid : ndarray
          The recursive ols residuals.
      rparams : ndarray
          The recursive ols parameter estimates.
      rypred : ndarray
          The recursive prediction of endogenous variable.
      rresid_standardized : ndarray
          The recursive residuals standardized so that N(0,sigma2) distributed,
          where sigma2 is the error variance.
      rresid_scaled : ndarray
          The recursive residuals normalize so that N(0,1) distributed.
      rcusum : ndarray
          The cumulative residuals for cusum test.
      rcusumci : ndarray
          The confidence interval for cusum test using a size of alpha.
  
      Notes
      -----
      It produces same recursive residuals as other version. This version updates
      the inverse of the X'X matrix and does not require matrix inversion during
      updating. looks efficient but no timing
  
      Confidence interval in Greene and Brown, Durbin and Evans is the same as
      in Ploberger after a little bit of algebra.
  
      References
      ----------
      jplv to check formulas, follows Harvey
      BigJudge 5.5.2b for formula for inverse(X'X) updating
      Greene section 7.5.2
  
      Brown, R. L., J. Durbin, and J. M. Evans. “Techniques for Testing the
      Constancy of Regression Relationships over Time.”
      Journal of the Royal Statistical Society. Series B (Methodological) 37,
      no. 2 (1975): 149-192.
      """
      if not isinstance(res, RegressionResultsWrapper):
          raise TypeError("res a regression results instance")
      y = res.model.endog
      x = res.model.exog
      order_by = array_like(
          order_by,
          "order_by",
          dtype="int",
          optional=True,
          ndim=1,
          shape=(y.shape[0],),
      )
      # intialize with skip observations
      if order_by is not None:
          x = x[order_by]
          y = y[order_by]
  
      nobs, nvars = x.shape
      if skip is None:
          skip = nvars
      rparams = np.nan * np.zeros((nobs, nvars))
      rresid = np.nan * np.zeros(nobs)
      rypred = np.nan * np.zeros(nobs)
      rvarraw = np.nan * np.zeros(nobs)
  
      x0 = x[:skip]
      if np.linalg.matrix_rank(x0) < x0.shape[1]:
          err_msg = """\
  "The initial regressor matrix, x[:skip], issingular. You must use a value of
  skip large enough to ensure that the first OLS estimator is well-defined.
  """
          raise ValueError(err_msg)
      y0 = y[:skip]
      # add Ridge to start (not in jplv)
      xtxi = np.linalg.inv(np.dot(x0.T, x0) + lamda * np.eye(nvars))
      xty = np.dot(x0.T, y0)  # xi * y   #np.dot(xi, y)
      beta = np.dot(xtxi, xty)
      rparams[skip - 1] = beta
      yipred = np.dot(x[skip - 1], beta)
      rypred[skip - 1] = yipred
      rresid[skip - 1] = y[skip - 1] - yipred
      rvarraw[skip - 1] = 1 + np.dot(
          x[skip - 1], np.dot(xtxi, x[skip - 1])
      )
      for i in range(skip, nobs):
          xi = x[i : i + 1, :]
          yi = y[i]
  
          # get prediction error with previous beta
          yipred = np.dot(xi, beta)
          rypred[i] = np.squeeze(yipred)
          residi = yi - yipred
          rresid[i] = np.squeeze(residi)
  
          # update beta and inverse(X'X)
          tmp = np.dot(xtxi, xi.T)
          ft = 1 + np.dot(xi, tmp)
  
          xtxi = xtxi - np.dot(tmp, tmp.T) / ft  # BigJudge equ 5.5.15
  
          beta = (
              beta + (tmp * residi / ft).ravel()
          )  # BigJudge equ 5.5.14
          rparams[i] = beta
          rvarraw[i] = np.squeeze(ft)
  
      rresid_scaled = rresid / np.sqrt(
          rvarraw
      )  # N(0,sigma2) distributed
      nrr = nobs - skip
      # sigma2 = rresid_scaled[skip-1:].var(ddof=1)  #var or sum of squares ?
      # Greene has var, jplv and Ploberger have sum of squares (Ass.:mean=0)
      # Gretl uses: by reverse engineering matching their numbers
      sigma2 = rresid_scaled[skip:].var(ddof=1)
      rresid_standardized = rresid_scaled / np.sqrt(
          sigma2
      )  # N(0,1) distributed
      rcusum = rresid_standardized[skip - 1 :].cumsum()
      # confidence interval points in Greene p136 looks strange. Cleared up
      # this assumes sum of independent standard normal, which does not take into
      # account that we make many tests at the same time
      if alpha == 0.90:
          a = 0.850
      elif alpha == 0.95:
          a = 0.948
      elif alpha == 0.99:
          a = 1.143
      else:
          raise ValueError("alpha can only be 0.9, 0.95 or 0.99")
  
      # following taken from Ploberger,
      # crit = a * np.sqrt(nrr)
      rcusumci = (
          a * np.sqrt(nrr)
          + 2 * a * np.arange(0, nobs - skip) / np.sqrt(nrr)
      ) * np.array([[-1.0], [+1.0]])
      return (
          rresid,
          rparams,
          rypred,
          rresid_standardized,
          rresid_scaled,
          rcusum,
          rcusumci,
      )
  
imports: from statsmodels.compat.pandas import deprecate_kwarg
  from collections.abc import Iterable
  import numpy
  import pandas
  from scipy import stats
  from statsmodels.regression.linear_model import OLS, RegressionResultsWrapper
  from statsmodels.stats._adnorm import anderson_statistic, normal_ad
  from statsmodels.stats._lilliefors import kstest_exponential, kstest_fit, kstest_normal, lilliefors
  from statsmodels.tools.validation import array_like, bool_like, dict_like, float_like, int_like, string_like
  from statsmodels.tsa.tsatools import lagmat
label_desc: Use image node_id 19 for calling a global function with example usage: recursive_olsresiduals(res, skip, lamda, alpha, order_by) and returns: rresid,  rparams,  rypred,  rresid_standardized,  rresid_scaled,  rcusum,  rcusumci