The R Project SVN R

/*
 *  R : A Computer Language for Statistical Data Analysis
 *  Copyright (C) 2002-2016   The R Core Team.
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, a copy is available at
 *  https://www.R-project.org/Licenses/
 */

#ifdef HAVE_CONFIG_H
# include <config.h>
#endif

#include <stdlib.h> // for abs
#include <string.h>

#include <R.h>
#include "ts.h"
#include "statsR.h" // for getListElement

#ifndef max
#define max(a,b) ((a < b)?(b):(a))
#endif
#ifndef min
#define min(a,b) ((a < b)?(a):(b))
#endif


/* 
  KalmanLike, internal to StructTS:
  .Call(C_KalmanLike, y, mod$Z, mod$a, mod$P, mod$T, mod$V, mod$h, mod$Pn,
        nit, FALSE, update)
  KalmanRun:
  .Call(C_KalmanLike, y, mod$Z, mod$a, mod$P, mod$T, mod$V, mod$h, mod$Pn,
        nit, TRUE, update)
*/

/* y vector length n of observations
   Z vector length p for observation equation y_t = Za_t +  eps_t
   a vector length p of initial state
   P p x p matrix for initial state uncertainty (contemparaneous)
   T  p x p transition matrix
   V  p x p = RQR'
   h = var(eps_t)
   anew used for a[t|t-1]
   Pnew used for P[t|t -1]
   M used for M = P[t|t -1]Z

   op is FALSE for KalmanLike, TRUE for KalmanRun.  
   The latter computes residuals and states and has 
   a more elaborate return value.

   Almost no checking here!
 */

SEXP
KalmanLike(SEXP sy, SEXP mod, SEXP sUP, SEXP op, SEXP update)
{
    int lop = asLogical(op);
    mod = PROTECT(duplicate(mod));

    SEXP sZ = getListElement(mod, "Z"), sa = getListElement(mod, "a"), 
    sP = getListElement(mod, "P"), sT = getListElement(mod, "T"), 
    sV = getListElement(mod, "V"), sh = getListElement(mod, "h"),
    sPn = getListElement(mod, "Pn");

    if (TYPEOF(sy) != REALSXP || TYPEOF(sZ) != REALSXP ||
    TYPEOF(sa) != REALSXP || TYPEOF(sP) != REALSXP ||
    TYPEOF(sPn) != REALSXP ||
    TYPEOF(sT) != REALSXP || TYPEOF(sV) != REALSXP)
    error(_("invalid argument type"));

    int n = LENGTH(sy), p = LENGTH(sa);
    double *y = REAL(sy), *Z = REAL(sZ), *T = REAL(sT), *V = REAL(sV),
    *P = REAL(sP), *a = REAL(sa), *Pnew = REAL(sPn), h = asReal(sh);

    double *anew = (double *) R_alloc(p, sizeof(double));
    double *M = (double *) R_alloc(p, sizeof(double));
    double *mm = (double *) R_alloc(p * p, sizeof(double));
    // These are only used if(lop), but avoid -Wall trouble
    SEXP ans = R_NilValue, resid = R_NilValue, states = R_NilValue;
    if(lop) {
    PROTECT(ans = allocVector(VECSXP, 3));
    SET_VECTOR_ELT(ans, 1, resid = allocVector(REALSXP, n));
    SET_VECTOR_ELT(ans, 2, states = allocMatrix(REALSXP, n, p));
    SEXP nm = PROTECT(allocVector(STRSXP, 3));
    SET_STRING_ELT(nm, 0, mkChar("values"));
    SET_STRING_ELT(nm, 1, mkChar("resid"));
    SET_STRING_ELT(nm, 2, mkChar("states"));
    setAttrib(ans, R_NamesSymbol, nm);
    UNPROTECT(1);
    }

    double sumlog = 0.0, ssq = 0.0;
    int nu = 0;
    for (int l = 0; l < n; l++) {
    for (int i = 0; i < p; i++) {
        double tmp = 0.0;
        for (int k = 0; k < p; k++)
        tmp += T[i + p * k] * a[k];
        anew[i] = tmp;
    }
    if (l > asInteger(sUP)) {
        for (int i = 0; i < p; i++)
        for (int j = 0; j < p; j++) {
            double tmp = 0.0;
            for (int k = 0; k < p; k++)
            tmp += T[i + p * k] * P[k + p * j];
            mm[i + p * j] = tmp;
        }
        for (int i = 0; i < p; i++)
        for (int j = 0; j < p; j++) {
            double tmp = V[i + p * j];
            for (int k = 0; k < p; k++)
            tmp += mm[i + p * k] * T[j + p * k];
            Pnew[i + p * j] = tmp;
        }
    }
    if (!ISNAN(y[l])) {
        nu++;
        double *rr = NULL /* -Wall */;
        if(lop) rr = REAL(resid);
        double resid0 = y[l];
        for (int i = 0; i < p; i++)
        resid0 -= Z[i] * anew[i];
        double gain = h;
        for (int i = 0; i < p; i++) {
        double tmp = 0.0;
        for (int j = 0; j < p; j++)
            tmp += Pnew[i + j * p] * Z[j];
        M[i] = tmp;
        gain += Z[i] * M[i];
        }
        ssq += resid0 * resid0 / gain;
        if(lop) rr[l] = resid0 / sqrt(gain);
        sumlog += log(gain);
        for (int i = 0; i < p; i++)
        a[i] = anew[i] + M[i] * resid0 / gain;
        for (int i = 0; i < p; i++)
        for (int j = 0; j < p; j++)
            P[i + j * p] = Pnew[i + j * p] - M[i] * M[j] / gain;
    } else {
        double *rr = NULL /* -Wall */;
        if(lop) rr = REAL(resid);
        for (int i = 0; i < p; i++)
        a[i] = anew[i];
        for (int i = 0; i < p * p; i++)
        P[i] = Pnew[i];
        if(lop) rr[l] = NA_REAL;
    }
    if(lop) {
        double *rs = REAL(states);
        for (int j = 0; j < p; j++) rs[l + n*j] = a[j];
    }
    }

    SEXP res = PROTECT(allocVector(REALSXP, 2));
    REAL(res)[0] = ssq/nu; REAL(res)[1] = sumlog/nu;
    if(lop) {
    SET_VECTOR_ELT(ans, 0, res);
    if(asLogical(update)) setAttrib(ans, install("mod"), mod);
    UNPROTECT(3);
    return ans;
    } else {
    if(asLogical(update)) setAttrib(res, install("mod"), mod);
    UNPROTECT(2);
    return res;
    }
}

SEXP
KalmanSmooth(SEXP sy, SEXP mod, SEXP sUP)
{
    SEXP sZ = getListElement(mod, "Z"), sa = getListElement(mod, "a"), 
    sP = getListElement(mod, "P"), sT = getListElement(mod, "T"), 
    sV = getListElement(mod, "V"), sh = getListElement(mod, "h"),
    sPn = getListElement(mod, "Pn");

    if (TYPEOF(sy) != REALSXP || TYPEOF(sZ) != REALSXP ||
    TYPEOF(sa) != REALSXP || TYPEOF(sP) != REALSXP ||
    TYPEOF(sT) != REALSXP || TYPEOF(sV) != REALSXP)
    error(_("invalid argument type"));

    SEXP ssa, ssP, ssPn, res, states = R_NilValue, sN;
    int n = LENGTH(sy), p = LENGTH(sa);
    double *y = REAL(sy), *Z = REAL(sZ), *a, *P,
    *T = REAL(sT), *V = REAL(sV), h = asReal(sh), *Pnew;
    double *at, *rt, *Pt, *gains, *resids, *Mt, *L, gn, *Nt;
    Rboolean var = TRUE;

    PROTECT(ssa = duplicate(sa)); a = REAL(ssa);
    PROTECT(ssP = duplicate(sP)); P = REAL(ssP);
    PROTECT(ssPn = duplicate(sPn)); Pnew = REAL(ssPn);

    PROTECT(res = allocVector(VECSXP, 2));
    SEXP nm = PROTECT(allocVector(STRSXP, 2));
    SET_STRING_ELT(nm, 0, mkChar("smooth"));
    SET_STRING_ELT(nm, 1, mkChar("var"));
    setAttrib(res, R_NamesSymbol, nm);
    UNPROTECT(1);
    SET_VECTOR_ELT(res, 0, states = allocMatrix(REALSXP, n, p));
    at = REAL(states);
    SET_VECTOR_ELT(res, 1, sN = allocVector(REALSXP, n*p*p));
    Nt = REAL(sN);

    double *anew, *mm, *M;
    anew = (double *) R_alloc(p, sizeof(double));
    M = (double *) R_alloc(p, sizeof(double));
    mm = (double *) R_alloc(p * p, sizeof(double));

    Pt = (double *) R_alloc(n * p * p, sizeof(double));
    gains = (double *) R_alloc(n, sizeof(double));
    resids = (double *) R_alloc(n, sizeof(double));
    Mt = (double *) R_alloc(n * p, sizeof(double));
    L = (double *) R_alloc(p * p, sizeof(double));

    for (int l = 0; l < n; l++) {
    for (int i = 0; i < p; i++) {
        double tmp = 0.0;
        for (int k = 0; k < p; k++)
        tmp += T[i + p * k] * a[k];
        anew[i] = tmp;
    }
    if (l > asInteger(sUP)) {
        for (int i = 0; i < p; i++)
        for (int j = 0; j < p; j++) {
            double tmp = 0.0;
            for (int k = 0; k < p; k++)
            tmp += T[i + p * k] * P[k + p * j];
            mm[i + p * j] = tmp;
        }
        for (int i = 0; i < p; i++)
        for (int j = 0; j < p; j++) {
            double tmp = V[i + p * j];
            for (int k = 0; k < p; k++)
            tmp += mm[i + p * k] * T[j + p * k];
            Pnew[i + p * j] = tmp;
        }
    }
    for (int i = 0; i < p; i++) at[l + n*i] = anew[i];
    for (int i = 0; i < p*p; i++) Pt[l + n*i] = Pnew[i];
    if (!ISNAN(y[l])) {
        double resid0 = y[l];
        for (int i = 0; i < p; i++)
        resid0 -= Z[i] * anew[i];
        double gain = h;
        for (int i = 0; i < p; i++) {
        double tmp = 0.0;
        for (int j = 0; j < p; j++)
            tmp += Pnew[i + j * p] * Z[j];
        Mt[l + n*i] = M[i] = tmp;
        gain += Z[i] * M[i];
        }
        gains[l] = gain;
        resids[l] = resid0;
        for (int i = 0; i < p; i++)
        a[i] = anew[i] + M[i] * resid0 / gain;
        for (int i = 0; i < p; i++)
        for (int j = 0; j < p; j++)
            P[i + j * p] = Pnew[i + j * p] - M[i] * M[j] / gain;
    } else {
        for (int i = 0; i < p; i++) {
        a[i] = anew[i];
        Mt[l + n * i] = 0.0;
        }
        for (int i = 0; i < p * p; i++)
        P[i] = Pnew[i];
        gains[l] = NA_REAL;
        resids[l] = NA_REAL;
    }
    }

    /* rt stores r_{t-1} */
    rt = (double *) R_alloc(n * p, sizeof(double));
    for (int l = n - 1; l >= 0; l--) {
    if (!ISNAN(gains[l])) {
        gn = 1/gains[l];
        for (int i = 0; i < p; i++)
        rt[l + n * i] = Z[i] * resids[l] * gn;
    } else {
        for (int i = 0; i < p; i++) rt[l + n * i] = 0.0;
        gn = 0.0;
    }

    if (var) {
        for (int i = 0; i < p; i++)
        for (int j = 0; j < p; j++)
            Nt[l + n*i + n*p*j] = Z[i] * Z[j] * gn;
    }

    if (l < n - 1) {
        /* compute r_{t-1} */
        for (int i = 0; i < p; i++)
        for (int j = 0; j < p; j++)
            mm[i + p * j] = ((i==j) ? 1:0) - Mt[l + n * i] * Z[j] * gn;
        for (int i = 0; i < p; i++)
        for (int j = 0; j < p; j++) {
            double tmp = 0.0;
            for (int k = 0; k < p; k++)
            tmp += T[i + p * k] * mm[k + p * j];
            L[i + p * j] = tmp;
        }
        for (int i = 0; i < p; i++) {
        double tmp = 0.0;
        for (int j = 0; j < p; j++)
            tmp += L[j + p * i] * rt[l + 1 + n * j];
        rt[l + n * i] += tmp;
        }
        if(var) { /* compute N_{t-1} */
        for (int i = 0; i < p; i++)
            for (int j = 0; j < p; j++) {
            double tmp = 0.0;
            for (int k = 0; k < p; k++)
                tmp += L[k + p * i] * Nt[l + 1 + n*k + n*p*j];
            mm[i + p * j] = tmp;
            }
        for (int i = 0; i < p; i++)
            for (int j = 0; j < p; j++) {
            double tmp = 0.0;
            for (int k = 0; k < p; k++)
                tmp += mm[i + p * k] * L[k + p * j];
            Nt[l + n*i + n*p*j] += tmp;
            }
        }
    }

    for (int i = 0; i < p; i++) {
        double tmp = 0.0;
        for (int j = 0; j < p; j++)
        tmp += Pt[l + n*i + n*p*j] * rt[l + n * j];
        at[l + n*i] += tmp;
    }
    }
    if (var)
    for (int l = 0; l < n; l++) {
        for (int i = 0; i < p; i++)
        for (int j = 0; j < p; j++) {
            double tmp = 0.0;
            for (int k = 0; k < p; k++)
            tmp += Pt[l + n*i + n*p*k] * Nt[l + n*k + n*p*j];
            mm[i + p * j] = tmp;
        }
        for (int i = 0; i < p; i++)
        for (int j = 0; j < p; j++) {
            double tmp = Pt[l + n*i + n*p*j];
            for (int k = 0; k < p; k++)
            tmp -= mm[i + p * k] * Pt[l + n*k + n*p*j];
            Nt[l + n*i + n*p*j] = tmp;
        }
    }
    UNPROTECT(4);
    return res;
}


SEXP
KalmanFore(SEXP nahead, SEXP mod, SEXP update)
{
    mod = PROTECT(duplicate(mod));
    SEXP sZ = getListElement(mod, "Z"), sa = getListElement(mod, "a"), 
    sP = getListElement(mod, "P"), sT = getListElement(mod, "T"), 
    sV = getListElement(mod, "V"), sh = getListElement(mod, "h");

    if (TYPEOF(sZ) != REALSXP ||
    TYPEOF(sa) != REALSXP || TYPEOF(sP) != REALSXP ||
    TYPEOF(sT) != REALSXP || TYPEOF(sV) != REALSXP)
    error(_("invalid argument type"));

    int  n = asInteger(nahead), p = LENGTH(sa);
    double *Z = REAL(sZ), *a = REAL(sa), *P = REAL(sP), *T = REAL(sT),
    *V = REAL(sV), h = asReal(sh);
    double *mm, *anew, *Pnew;

    anew = (double *) R_alloc(p, sizeof(double));
    Pnew = (double *) R_alloc(p * p, sizeof(double));
    mm = (double *) R_alloc(p * p, sizeof(double));
    SEXP res, forecasts, se;
    PROTECT(res = allocVector(VECSXP, 2));
    SET_VECTOR_ELT(res, 0, forecasts = allocVector(REALSXP, n));
    SET_VECTOR_ELT(res, 1, se = allocVector(REALSXP, n));
    {
    SEXP nm = PROTECT(allocVector(STRSXP, 2));
    SET_STRING_ELT(nm, 0, mkChar("pred"));
    SET_STRING_ELT(nm, 1, mkChar("var"));
    setAttrib(res, R_NamesSymbol, nm);
    UNPROTECT(1);
    }
    for (int l = 0; l < n; l++) {
    double fc = 0.0;
    for (int i = 0; i < p; i++) {
        double tmp = 0.0;
        for (int k = 0; k < p; k++)
        tmp += T[i + p * k] * a[k];
        anew[i] = tmp;
        fc += tmp * Z[i];
    }
    for (int i = 0; i < p; i++)
        a[i] = anew[i];
    REAL(forecasts)[l] = fc;

    for (int i = 0; i < p; i++)
        for (int j = 0; j < p; j++) {
        double tmp = 0.0;
        for (int k = 0; k < p; k++)
            tmp += T[i + p * k] * P[k + p * j];
        mm[i + p * j] = tmp;
        }
    for (int i = 0; i < p; i++)
        for (int j = 0; j < p; j++) {
        double tmp = V[i + p * j];
        for (int k = 0; k < p; k++)
            tmp += mm[i + p * k] * T[j + p * k];
        Pnew[i + p * j] = tmp;
        }
    double tmp = h;
    for (int i = 0; i < p; i++)
        for (int j = 0; j < p; j++) {
        P[i + j * p] = Pnew[i + j * p];
        tmp += Z[i] * Z[j] * P[i + j * p];
        }
    REAL(se)[l] = tmp;
    }
    if(asLogical(update)) setAttrib(res, install("mod"), mod);
    UNPROTECT(2);
    return res;
}


static void partrans(int p, double *raw, double *new)
{
    int j, k;
    double a, work[100];

    if(p > 100) error(_("can only transform 100 pars in arima0"));

    /* Step one: map (-Inf, Inf) to (-1, 1) via tanh
       The parameters are now the pacf phi_{kk} */
    for(j = 0; j < p; j++) work[j] = new[j] = tanh(raw[j]);
    /* Step two: run the Durbin-Levinson recursions to find phi_{j.},
       j = 2, ..., p and phi_{p.} are the autoregression coefficients */
    for(j = 1; j < p; j++) {
    a = new[j];
    for(k = 0; k < j; k++)
        work[k] -= a * new[j - k - 1];
    for(k = 0; k < j; k++) new[k] = work[k];
    }
}

SEXP ARIMA_undoPars(SEXP sin, SEXP sarma)
{
    int *arma = INTEGER(sarma), mp = arma[0], mq = arma[1], msp = arma[2],
    v, n = LENGTH(sin);
    double *params, *in = REAL(sin);
    SEXP res = allocVector(REALSXP, n);

    params = REAL(res);
    for (int i = 0; i < n; i++) params[i] = in[i];
    if (mp > 0) partrans(mp, in, params);
    v = mp + mq;
    if (msp > 0) partrans(msp, in + v, params + v);
    return res;
}


SEXP ARIMA_transPars(SEXP sin, SEXP sarma, SEXP strans)
{
    int *arma = INTEGER(sarma), trans = asLogical(strans);
    int mp = arma[0], mq = arma[1], msp = arma[2], msq = arma[3],
    ns = arma[4], i, j, p = mp + ns * msp, q = mq + ns * msq, v;
    double *in = REAL(sin), *params = REAL(sin), *phi, *theta;
    SEXP res, sPhi, sTheta;

    PROTECT(res = allocVector(VECSXP, 2));
    SET_VECTOR_ELT(res, 0, sPhi = allocVector(REALSXP, p));
    SET_VECTOR_ELT(res, 1, sTheta = allocVector(REALSXP, q));
    phi = REAL(sPhi);
    theta = REAL(sTheta);
    if (trans) {
    int n = mp + mq + msp + msq;

    params = (double *) R_alloc(n, sizeof(double));
    for (i = 0; i < n; i++) params[i] = in[i];
    if (mp > 0) partrans(mp, in, params);
    v = mp + mq;
    if (msp > 0) partrans(msp, in + v, params + v);
    }
    if (ns > 0) {
    /* expand out seasonal ARMA models */
    for (i = 0; i < mp; i++) phi[i] = params[i];
    for (i = 0; i < mq; i++) theta[i] = params[i + mp];
    for (i = mp; i < p; i++) phi[i] = 0.0;
    for (i = mq; i < q; i++) theta[i] = 0.0;
    for (j = 0; j < msp; j++) {
        phi[(j + 1) * ns - 1] += params[j + mp + mq];
        for (i = 0; i < mp; i++)
        phi[(j + 1) * ns + i] -= params[i] * params[j + mp + mq];
    }
    for (j = 0; j < msq; j++) {
        theta[(j + 1) * ns - 1] += params[j + mp + mq + msp];
        for (i = 0; i < mq; i++)
        theta[(j + 1) * ns + i] += params[i + mp] *
            params[j + mp + mq + msp];
    }
    } else {
    for (i = 0; i < mp; i++) phi[i] = params[i];
    for (i = 0; i < mq; i++) theta[i] = params[i + mp];
    }
    UNPROTECT(1);
    return res;
}

#if !defined(atanh) && defined(HAVE_DECL_ATANH) && !HAVE_DECL_ATANH
extern double atanh(double x);
#endif
static void invpartrans(int p, double *phi, double *new)
{
    int j, k;
    double a, work[100];

    if(p > 100) error(_("can only transform 100 pars in arima0"));

    for(j = 0; j < p; j++) work[j] = new[j] = phi[j];
    /* Run the Durbin-Levinson recursions backwards
       to find the PACF phi_{j.} from the autoregression coefficients */
    for(j = p - 1; j > 0; j--) {
    a = new[j];
    for(k = 0; k < j; k++)
        work[k]  = (new[k] + a * new[j - k - 1]) / (1 - a * a);
    for(k = 0; k < j; k++) new[k] = work[k];
    }
    for(j = 0; j < p; j++) new[j] = atanh(new[j]);
}

SEXP ARIMA_Invtrans(SEXP in, SEXP sarma)
{
    int *arma = INTEGER(sarma), mp = arma[0], mq = arma[1], msp = arma[2],
    i, v, n = LENGTH(in);
    SEXP y = allocVector(REALSXP, n);
    double *raw = REAL(in), *new = REAL(y);

    for(i = 0; i < n; i++) new[i] = raw[i];
    if (mp > 0) invpartrans(mp, raw, new);
    v = mp + mq;
    if (msp > 0) invpartrans(msp, raw + v, new + v);
    return y;
}

#define eps 1e-3
SEXP ARIMA_Gradtrans(SEXP in, SEXP sarma)
{
    int *arma = INTEGER(sarma), mp = arma[0], mq = arma[1], msp = arma[2],
    n = LENGTH(in);
    SEXP y = allocMatrix(REALSXP, n, n);
    double *raw = REAL(in), *A = REAL(y), w1[100], w2[100], w3[100];

    for (int i = 0; i < n; i++)
    for (int j = 0; j < n; j++)
        A[i + j*n] = (i == j);
    if(mp > 0) {
    for (int i = 0; i < mp; i++) w1[i] = raw[i];
    partrans(mp, w1, w2);
    for (int i = 0; i < mp; i++) {
        w1[i] += eps;
        partrans(mp, w1, w3);
        for (int j = 0; j < mp; j++) A[i + j*n] = (w3[j] - w2[j])/eps;
        w1[i] -= eps;
    }
    }
    if(msp > 0) {
    int v = mp + mq;
    for (int i = 0; i < msp; i++) w1[i] = raw[i + v];
    partrans(msp, w1, w2);
    for(int i = 0; i < msp; i++) {
        w1[i] += eps;
        partrans(msp, w1, w3);
        for(int j = 0; j < msp; j++)
        A[i + v + (j+v)*n] = (w3[j] - w2[j])/eps;
        w1[i] -= eps;
    }
    }
    return y;
}


SEXP
ARIMA_Like(SEXP sy, SEXP mod, SEXP sUP, SEXP giveResid)
{
    SEXP sPhi = getListElement(mod, "phi"), 
    sTheta = getListElement(mod, "theta"), 
    sDelta = getListElement(mod, "Delta"),
    sa = getListElement(mod, "a"),
    sP = getListElement(mod, "P"),
    sPn = getListElement(mod, "Pn");

    if (TYPEOF(sPhi) != REALSXP || TYPEOF(sTheta) != REALSXP ||
    TYPEOF(sDelta) != REALSXP || TYPEOF(sa) != REALSXP ||
    TYPEOF(sP) != REALSXP || TYPEOF(sPn) != REALSXP)
    error(_("invalid argument type"));

    SEXP res, nres, sResid = R_NilValue;
    int n = LENGTH(sy), rd = LENGTH(sa), p = LENGTH(sPhi),
    q = LENGTH(sTheta), d = LENGTH(sDelta), r = rd - d;
    double *y = REAL(sy), *a = REAL(sa), *P = REAL(sP), *Pnew = REAL(sPn);
    double *phi = REAL(sPhi), *theta = REAL(sTheta), *delta = REAL(sDelta);
    double sumlog = 0.0, ssq = 0, *anew, *mm = NULL, *M;
    int nu = 0;
    Rboolean useResid = asLogical(giveResid);
    double *rsResid = NULL /* -Wall */;

    anew = (double *) R_alloc(rd, sizeof(double));
    M = (double *) R_alloc(rd, sizeof(double));
    if (d > 0) mm = (double *) R_alloc(rd * rd, sizeof(double));

    if (useResid) {
    PROTECT(sResid = allocVector(REALSXP, n));
    rsResid = REAL(sResid);
    }

    for (int l = 0; l < n; l++) {
    for (int i = 0; i < r; i++) {
        double tmp = (i < r - 1) ? a[i + 1] : 0.0;
        if (i < p) tmp += phi[i] * a[0];
        anew[i] = tmp;
    }
    if (d > 0) {
        for (int i = r + 1; i < rd; i++) anew[i] = a[i - 1];
        double tmp = a[0];
        for (int i = 0; i < d; i++) tmp += delta[i] * a[r + i];
        anew[r] = tmp;
    }
    if (l > asInteger(sUP)) {
        if (d == 0) {
        for (int i = 0; i < r; i++) {
            double vi = 0.0;
            if (i == 0) vi = 1.0; else if (i - 1 < q) vi = theta[i - 1];
            for (int j = 0; j < r; j++) {
            double tmp = 0.0;
            if (j == 0) tmp = vi; else if (j - 1 < q) tmp = vi * theta[j - 1];
            if (i < p && j < p) tmp += phi[i] * phi[j] * P[0];
            if (i < r - 1 && j < r - 1) tmp += P[i + 1 + r * (j + 1)];
            if (i < p && j < r - 1) tmp += phi[i] * P[j + 1];
            if (j < p && i < r - 1) tmp += phi[j] * P[i + 1];
            Pnew[i + r * j] = tmp;
            }
        }
        } else {
        /* mm = TP */
        for (int i = 0; i < r; i++)
            for (int j = 0; j < rd; j++) {
            double tmp = 0.0;
            if (i < p) tmp += phi[i] * P[rd * j];
            if (i < r - 1) tmp += P[i + 1 + rd * j];
            mm[i + rd * j] = tmp;
            }
        for (int j = 0; j < rd; j++) {
            double tmp = P[rd * j];
            for (int k = 0; k < d; k++)
            tmp += delta[k] * P[r + k + rd * j];
            mm[r + rd * j] = tmp;
        }
        for (int i = 1; i < d; i++)
            for (int j = 0; j < rd; j++)
            mm[r + i + rd * j] = P[r + i - 1 + rd * j];

        /* Pnew = mmT' */
        for (int i = 0; i < r; i++)
            for (int j = 0; j < rd; j++) {
            double tmp = 0.0;
            if (i < p) tmp += phi[i] * mm[j];
            if (i < r - 1) tmp += mm[rd * (i + 1) + j];
            Pnew[j + rd * i] = tmp;
            }
        for (int j = 0; j < rd; j++) {
            double tmp = mm[j];
            for (int k = 0; k < d; k++)
            tmp += delta[k] * mm[rd * (r + k) + j];
            Pnew[rd * r + j] = tmp;
        }
        for (int i = 1; i < d; i++)
            for (int j = 0; j < rd; j++)
            Pnew[rd * (r + i) + j] = mm[rd * (r + i - 1) + j];
        /* Pnew <- Pnew + (1 theta) %o% (1 theta) */
        for (int i = 0; i <= q; i++) {
            double vi = (i == 0) ? 1. : theta[i - 1];
            for (int j = 0; j <= q; j++)
            Pnew[i + rd * j] += vi * ((j == 0) ? 1. : theta[j - 1]);
        }
        }
    }
    if (!ISNAN(y[l])) {
        double resid = y[l] - anew[0];
        for (int i = 0; i < d; i++)
        resid -= delta[i] * anew[r + i];

        for (int i = 0; i < rd; i++) {
        double tmp = Pnew[i];
        for (int j = 0; j < d; j++)
            tmp += Pnew[i + (r + j) * rd] * delta[j];
        M[i] = tmp;
        }

        double gain = M[0];
        for (int j = 0; j < d; j++) gain += delta[j] * M[r + j];
        if(gain < 1e4) {
        nu++;
        ssq += resid * resid / gain;
        sumlog += log(gain);
        }
        if (useResid) rsResid[l] = resid / sqrt(gain);
        for (int i = 0; i < rd; i++)
        a[i] = anew[i] + M[i] * resid / gain;
        for (int i = 0; i < rd; i++)
        for (int j = 0; j < rd; j++)
            P[i + j * rd] = Pnew[i + j * rd] - M[i] * M[j] / gain;
    } else {
        for (int i = 0; i < rd; i++) a[i] = anew[i];
        for (int i = 0; i < rd * rd; i++) P[i] = Pnew[i];
        if (useResid) rsResid[l] = NA_REAL;
    }
    }

    if (useResid) {
    PROTECT(res = allocVector(VECSXP, 3));
    SET_VECTOR_ELT(res, 0, nres = allocVector(REALSXP, 3));
    REAL(nres)[0] = ssq;
    REAL(nres)[1] = sumlog;
    REAL(nres)[2] = (double) nu;
    SET_VECTOR_ELT(res, 1, sResid);
    UNPROTECT(2);
    return res;
    } else {
    nres = allocVector(REALSXP, 3);
    REAL(nres)[0] = ssq;
    REAL(nres)[1] = sumlog;
    REAL(nres)[2] = (double) nu;
    return nres;
    }
}

/* do differencing here */
/* arma is p, q, sp, sq, ns, d, sd */
SEXP
ARIMA_CSS(SEXP sy, SEXP sarma, SEXP sPhi, SEXP sTheta,
      SEXP sncond, SEXP giveResid)
{
    SEXP res, sResid = R_NilValue;
    double ssq = 0.0, *y = REAL(sy), tmp;
    double *phi = REAL(sPhi), *theta = REAL(sTheta), *w, *resid;
    int n = LENGTH(sy), *arma = INTEGER(sarma), p = LENGTH(sPhi),
    q = LENGTH(sTheta), ncond = asInteger(sncond);
    int ns, nu = 0;
    Rboolean useResid = asLogical(giveResid);

    w = (double *) R_alloc(n, sizeof(double));
    for (int l = 0; l < n; l++) w[l] = y[l];
    for (int i = 0; i < arma[5]; i++)
    for (int l = n - 1; l > 0; l--) w[l] -= w[l - 1];
    ns = arma[4];
    for (int i = 0; i < arma[6]; i++)
    for (int l = n - 1; l >= ns; l--) w[l] -= w[l - ns];

    PROTECT(sResid = allocVector(REALSXP, n));
    resid = REAL(sResid);
    if (useResid) for (int l = 0; l < ncond; l++) resid[l] = 0;

    for (int l = ncond; l < n; l++) {
    tmp = w[l];
    for (int j = 0; j < p; j++) tmp -= phi[j] * w[l - j - 1];
    for (int j = 0; j < min(l - ncond, q); j++)
        tmp -= theta[j] * resid[l - j - 1];
    resid[l] = tmp;
    if (!ISNAN(tmp)) {
        nu++;
        ssq += tmp * tmp;
    }
    }
    if (useResid) {
    PROTECT(res = allocVector(VECSXP, 2));
    SET_VECTOR_ELT(res, 0, ScalarReal(ssq / (double) (nu)));
    SET_VECTOR_ELT(res, 1, sResid);
    UNPROTECT(2);
    return res;
    } else {
    UNPROTECT(1);
    return ScalarReal(ssq / (double) (nu));
    }
}

SEXP TSconv(SEXP a, SEXP b)
{
    int na, nb, nab;
    SEXP ab;
    double *ra, *rb, *rab;

    PROTECT(a = coerceVector(a, REALSXP));
    PROTECT(b = coerceVector(b, REALSXP));
    na = LENGTH(a);
    nb = LENGTH(b);
    nab = na + nb - 1;
    PROTECT(ab = allocVector(REALSXP, nab));
    ra = REAL(a); rb = REAL(b); rab = REAL(ab);
    for (int i = 0; i < nab; i++) rab[i] = 0.0;
    for (int i = 0; i < na; i++)
    for (int j = 0; j < nb; j++)
        rab[i + j] += ra[i] * rb[j];
    UNPROTECT(3);
    return (ab);
}

/* based on code from AS154 */

static void
inclu2(size_t np, double *xnext, double *xrow, double ynext,
       double *d, double *rbar, double *thetab)
{
    double cbar, sbar, di, xi, xk, rbthis, dpi;
    size_t i, k, ithisr;

/*   This subroutine updates d, rbar, thetab by the inclusion
     of xnext and ynext. */

    for (i = 0; i < np; i++) xrow[i] = xnext[i];

    for (ithisr = 0, i = 0; i < np; i++) {
    if (xrow[i] != 0.0) {
        xi = xrow[i];
        di = d[i];
        dpi = di + xi * xi;
        d[i] = dpi;
        cbar = di / dpi;
        sbar = xi / dpi;
        for (k = i + 1; k < np; k++) {
        xk = xrow[k];
        rbthis = rbar[ithisr];
        xrow[k] = xk - xi * rbthis;
        rbar[ithisr++] = cbar * rbthis + sbar * xk;
        }
        xk = ynext;
        ynext = xk - xi * thetab[i];
        thetab[i] = cbar * thetab[i] + sbar * xk;
        if (di == 0.0) return;
    } else
        ithisr = ithisr + np - i - 1;
    }
}

#ifdef DEBUG_Q0bis
# include <R_ext/Print.h>
  double chk_V(double v[], char* nm, int jj, int len) { 
    // len = length(<vector>)  <==> index must be in  {0, len-1}
    if(jj < 0 || jj >= len)
    REprintf(" %s[%2d]\n", nm, jj); 
    return(v[jj]); 
  }
#endif

/*
  Matwey V. Kornilov's implementation of algorithm by
  Dr. Raphael Rossignol
  See https://bugs.r-project.org/show_bug.cgi?id=14682 for details.
*/
SEXP getQ0bis(SEXP sPhi, SEXP sTheta, SEXP sTol)
{
    SEXP res;
    int p = LENGTH(sPhi), q = LENGTH(sTheta);
    double *phi = REAL(sPhi), *theta = REAL(sTheta); // tol = REAL(sTol)[0];

    int i,j, r = max(p, q + 1);

    /* Final result is block product 
     *   Q0 = A1 SX A1^T + A1 SXZ A2^T + (A1 SXZ A2^T)^T + A2 A2^T ,
     * where A1 [i,j] = phi[i+j],
     *       A2 [i,j] = ttheta[i+j],  and SX, SXZ are defined below */
    PROTECT(res = allocMatrix(REALSXP, r, r));
    double *P = REAL(res);

    /* Clean P */
    Memzero(P, r*r);

#ifdef DEBUG_Q0bis
#define _ttheta(j) chk_V(ttheta, "ttheta", j, q+1)// was  r
#define _tphi(j)   chk_V(tphi,   "tphi",   j, p+1)
#define _rrz(j)    chk_V(rrz,    "rrz",    j, q)
#else
#define _ttheta(j) ttheta[j]
#define _tphi(j) tphi[j]
#define _rrz(j)  rrz [j]
#endif

    double *ttheta = (double *) R_alloc(q + 1, sizeof(double));
    /* Init ttheta = c(1, theta) */
    ttheta[0] = 1.;
    for (i = 1; i < q + 1; ++i) ttheta[i] = theta[i - 1];

    if( p > 0 ) {
    int r2 = max(p + q, p + 1);
    SEXP sgam = PROTECT(allocMatrix(REALSXP, r2, r2)),
        sg = PROTECT(allocVector(REALSXP, r2));
    double *gam = REAL(sgam);
    double *g = REAL(sg);
    double *tphi = (double *) R_alloc(p + 1, sizeof(double));
    /* Init tphi = c(1, -phi) */
    tphi[0] = 1.;
    for (i = 1; i < p + 1; ++i) tphi[i] = -phi[i - 1];

    /* Compute the autocovariance function of U, the AR part of X */

    /* Gam := C1 + C2 ; initialize */
    Memzero(gam, r2*r2);

    /* C1[E] */
    for (j = 0; j < r2; ++j)
        for (i = j; i < r2 && i - j < p + 1; ++i)
        gam[j*r2 + i] += _tphi(i-j);

    /* C2[E] */
    for (i = 0; i < r2; ++i)
        for (j = 1; j < r2 && i + j < p + 1; ++j)
        gam[j*r2 + i] += _tphi(i+j);

    /* Initialize g = (1 0 0 .... 0) */
    g[0] = 1.;
    for (i = 1; i < r2; ++i)
        g[i] = 0.;

    /* rU = solve(Gam, g)  -> solve.default() -> .Internal(La_solve, .,)
     * --> fiddling with R-objects -> C and then F77_CALL(.) of dgesv, dlange, dgecon
     * FIXME: call these directly here, possibly even use 'info' instead of error(.)
     * e.g., in case of exact singularity.
     */
    SEXP callS = PROTECT(lang4(install("solve.default"), sgam, sg, sTol)),
        su = PROTECT(eval(callS, R_BaseEnv));
    double *u = REAL(su);
    /* SX = A SU A^T */
    /* A[i,j]  = ttheta[j-i] */
    /* SU[i,j] = u[abs(i-j)] */
    /* Q0 += ( A1 SX A1^T == A1 A SU A^T A1^T) */
    // (relying on good compiler optimization here:)
    for (i = 0; i < r; ++i)
        for (j = i; j < r; ++j)
        for (int k = 0; i + k < p; ++k)
            for (int L = k; L - k < q + 1; ++L)
            for (int m = 0; j + m < p; ++m)
                for (int n = m; n - m < q + 1; ++n)
                P[r*i + j] += phi[i + k] * phi[j + m] *
                    _ttheta(L - k) * _ttheta(n - m) * u[abs(L - n)];
    UNPROTECT(4);
    /* Compute correlation matrix between X and Z */
    /* forwardsolve(C1, g) */
    /* C[i,j] = tphi[i-j] */
    /* g[i] = _ttheta(i) */
    double *rrz = (double *) R_alloc(q, sizeof(double));
    if(q > 0) {
        for (i = 0; i < q; ++i) {
        rrz[i] = _ttheta(i);
        for (j = max(0, i - p); j < i; ++j)
            rrz[i] -= _rrz(j) * _tphi(i-j);
        }
    }

    /* Q0 += A1 SXZ A2^T + (A1 SXZ A2^T)^T */
    /* SXZ[i,j] = rrz[j-i-1], j > 0 */
    for (i = 0; i < r; ++i)
        for (j = i; j < r; ++j) {
        int k, L;
        for (k = 0; i + k < p; ++k)
            for (L = k+1; j + L < q + 1; ++L)
            P[r*i + j] += phi[i + k] * _ttheta(j + L) * _rrz(L - k - 1);
        for (k = 0; j + k < p; ++k)
            for (L = k+1; i + L < q + 1; ++L)
            P[r*i + j] += phi[j + k] * _ttheta(i + L) * _rrz(L - k - 1);
        }
    } // end if(p > 0)

    /* Q0 += A2 A2^T */
    for (i = 0; i < r; ++i)
    for (j = i; j < r; ++j)
        for (int k = 0; j + k < q + 1; ++k)
         P[r*i + j] += _ttheta(i + k) * _ttheta(j + k);

    /* Symmetrize result */
    for (i = 0; i < r; ++i)
    for (j = i+1; j < r; ++j)
        P[r*j + i] = P[r*i + j];

    UNPROTECT(1);
    return res;
}

SEXP getQ0(SEXP sPhi, SEXP sTheta)
{
    SEXP res;
    int  p = LENGTH(sPhi), q = LENGTH(sTheta);
    double *phi = REAL(sPhi), *theta = REAL(sTheta);

    /* thetab[np], xnext[np], xrow[np].  rbar[rbar] */
    /* NB: nrbar could overflow */
    int r = max(p, q + 1);
    size_t np = r * (r + 1) / 2, nrbar = np * (np - 1) / 2, npr, npr1;
    size_t indi, indj, indn, i, j, ithisr, ind, ind1, ind2, im, jm;


    /* This is the limit using an int index.  We could use
       size_t and get more on a 64-bit system,
       but there seems no practical need. */
    if(r > 350) error(_("maximum supported lag is 350"));
    double *xnext, *xrow, *rbar, *thetab, *V;
    xnext = (double *) R_alloc(np, sizeof(double));
    xrow = (double *) R_alloc(np, sizeof(double));
    rbar = (double *) R_alloc(nrbar, sizeof(double));
    thetab = (double *) R_alloc(np, sizeof(double));
    V = (double *) R_alloc(np, sizeof(double));
    for (ind = 0, j = 0; j < r; j++) {
    double vj = 0.0;
    if (j == 0) vj = 1.0; else if (j - 1 < q) vj = theta[j - 1];
    for (i = j; i < r; i++) {
        double vi = 0.0;
        if (i == 0) vi = 1.0; else if (i - 1 < q) vi = theta[i - 1];
        V[ind++] = vi * vj;
    }
    }

    PROTECT(res = allocMatrix(REALSXP, r, r));
    double *P = REAL(res);

    if (r == 1) {
    if (p == 0) P[0] = 1.0; // PR#16419
    else P[0] = 1.0 / (1.0 - phi[0] * phi[0]);
    UNPROTECT(1);
    return res;
    }
    if (p > 0) {
/*      The set of equations s * vec(P0) = vec(v) is solved for
    vec(P0).  s is generated row by row in the array xnext.  The
    order of elements in P is changed, so as to bring more leading
    zeros into the rows of s. */

    for (i = 0; i < nrbar; i++) rbar[i] = 0.0;
    for (i = 0; i < np; i++) {
        P[i] = 0.0;
        thetab[i] = 0.0;
        xnext[i] = 0.0;
    }
    ind = 0;
    ind1 = -1;
    npr = np - r;
    npr1 = npr + 1;
    indj = npr;
    ind2 = npr - 1;
    for (j = 0; j < r; j++) {
        double phij = (j < p) ? phi[j] : 0.0;
        xnext[indj++] = 0.0;
        indi = npr1 + j;
        for (i = j; i < r; i++) {
        double ynext = V[ind++];
        double phii = (i < p) ? phi[i] : 0.0;
        if (j != r - 1) {
            xnext[indj] = -phii;
            if (i != r - 1) {
            xnext[indi] -= phij;
            xnext[++ind1] = -1.0;
            }
        }
        xnext[npr] = -phii * phij;
        if (++ind2 >= np) ind2 = 0;
        xnext[ind2] += 1.0;
        inclu2(np, xnext, xrow, ynext, P, rbar, thetab);
        xnext[ind2] = 0.0;
        if (i != r - 1) {
            xnext[indi++] = 0.0;
            xnext[ind1] = 0.0;
        }
        }
    }

    ithisr = nrbar - 1;
    im = np - 1;
    for (i = 0; i < np; i++) {
        double bi = thetab[im];
        for (jm = np - 1, j = 0; j < i; j++)
        bi -= rbar[ithisr--] * P[jm--];
        P[im--] = bi;
    }

/*        now re-order p. */

    ind = npr;
    for (i = 0; i < r; i++) xnext[i] = P[ind++];
    ind = np - 1;
    ind1 = npr - 1;
    for (i = 0; i < npr; i++) P[ind--] = P[ind1--];
    for (i = 0; i < r; i++) P[i] = xnext[i];
    } else {

/* P0 is obtained by backsubstitution for a moving average process. */

    indn = np;
    ind = np;
    for (i = 0; i < r; i++)
        for (j = 0; j <= i; j++) {
        --ind;
        P[ind] = V[ind];
        if (j != 0) P[ind] += P[--indn];
        }
    }
    /* now unpack to a full matrix */
    for (i = r - 1, ind = np; i > 0; i--)
    for (j = r - 1; j >= i; j--)
        P[r * i + j] = P[--ind];
    for (i = 0; i < r - 1; i++)
    for (j = i + 1; j < r; j++)
        P[i + r * j] = P[j + r * i];
    UNPROTECT(1);
    return res;
}