R stats/arima.c

arima.c

/*
 *  R : A Computer Language for Statistical Data Analysis
 *  Copyright (C) 2002-2006   The R Development 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
 *  http://www.r-project.org/Licenses/
 */

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

#if !defined(atanh) && defined(HAVE_DECL_ATANH) && !HAVE_DECL_ATANH
extern double atanh(double x);
#endif

/* do this first to get the right options for math.h */
#include <R_ext/Arith.h>

#include <R.h>
#include "ts.h"

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


/* 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

   No checking here!
 */
SEXP
KalmanLike(SEXP sy, SEXP sZ, SEXP sa, SEXP sP, SEXP sT,
	   SEXP sV, SEXP sh, SEXP sPn, SEXP sUP, SEXP op, SEXP fast)
{
    SEXP res, ans = R_NilValue, resid = R_NilValue, states = R_NilValue;
    int n, p, lop = asLogical(op);
    double *y, *Z, *a, *P, *T, *V, h = asReal(sh), *Pnew;
    double sumlog = 0.0, ssq = 0, resid0, gain, tmp, *anew, *mm, *M;
    int i, j, k, l;

    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"));
    n = LENGTH(sy); p = LENGTH(sa);
    y = REAL(sy); Z = REAL(sZ); T = REAL(sT); V = REAL(sV);

    /* Avoid modifying arguments unless fast=TRUE */
    if (!LOGICAL(fast)[0]){
	    PROTECT(sP = duplicate(sP));
	    PROTECT(sa = duplicate(sa));
	    PROTECT(sPn = duplicate(sPn));
    }
    P = REAL(sP); a = REAL(sa); Pnew = REAL(sPn);

    anew = (double *) R_alloc(p, sizeof(double));
    M = (double *) R_alloc(p, sizeof(double));
    mm = (double *) R_alloc(p * p, sizeof(double));
    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));
    }
    for (l = 0; l < n; l++) {
	for (i = 0; i < p; i++) {
	    tmp = 0.0;
	    for (k = 0; k < p; k++)
		tmp += T[i + p * k] * a[k];
	    anew[i] = tmp;
	}
	if (l > asInteger(sUP)) {
	    for (i = 0; i < p; i++)
		for (j = 0; j < p; j++) {
		    tmp = 0.0;
		    for (k = 0; k < p; k++)
			tmp += T[i + p * k] * P[k + p * j];
		    mm[i + p * j] = tmp;
		}
	    for (i = 0; i < p; i++)
		for (j = 0; j < p; j++) {
		    tmp = V[i + p * j];
		    for (k = 0; k < p; k++)
			tmp += mm[i + p * k] * T[j + p * k];
		    Pnew[i + p * j] = tmp;
		}
	}
	if (!ISNAN(y[l])) {
	    double *rr = NULL /* -Wall */;
	    if(lop) rr = REAL(resid);
	    resid0 = y[l];
	    for (i = 0; i < p; i++)
		resid0 -= Z[i] * anew[i];
	    gain = h;
	    for (i = 0; i < p; i++) {
		tmp = 0.0;
		for (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 (i = 0; i < p; i++)
		a[i] = anew[i] + M[i] * resid0 / gain;
	    for (i = 0; i < p; i++)
		for (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 (i = 0; i < p; i++)
		a[i] = anew[i];
	    for (i = 0; i < p * p; i++)
		P[i] = Pnew[i];
	    if(lop) rr[l] = NA_REAL;
	}
	if(lop) {
	    double *rs = REAL(states);
	    for(j = 0; j < p; j++) rs[l + n*j] = a[j];
	}
    }

    if(lop) {
	SET_VECTOR_ELT(ans, 0, res=allocVector(REALSXP, 2));
	REAL(res)[0] = ssq; REAL(res)[1] = sumlog;
	UNPROTECT(1);
	if (!LOGICAL(fast)[0])
	    UNPROTECT(3);
	return ans;
    } else {
	res = allocVector(REALSXP, 2);
	REAL(res)[0] = ssq; REAL(res)[1] = sumlog;
	if (!LOGICAL(fast)[0])
	    UNPROTECT(3);
	return res;
    }
}

SEXP
KalmanSmooth(SEXP sy, SEXP sZ, SEXP sa, SEXP sP, SEXP sT,
	     SEXP sV, SEXP sh, SEXP sPn, SEXP sUP)
{
    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 resid0, gain, tmp, *anew, *mm, *M;
    double *at, *rt, *Pt, *gains, *resids, *Mt, *L, gn, *Nt;
    int i, j, k, l;
    Rboolean var = TRUE;

    /* It would be better to check types before using LENGTH and REAL
       on these, but should still work this way.  LT */
    if (TYPEOF(sy) != REALSXP || TYPEOF(sZ) != REALSXP ||
	TYPEOF(sa) != REALSXP || TYPEOF(sP) != REALSXP ||
	TYPEOF(sT) != REALSXP || TYPEOF(sV) != REALSXP)
	error(_("invalid argument type"));

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

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

	if (var) {
	    for (i = 0; i < p; i++)
		for (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 (i = 0; i < p; i++)
		for (j = 0; j < p; j++)
		    mm[i + p * j] = ((i==j)?1:0) - Mt[l + n * i] * Z[j] * gn;
	    for (i = 0; i < p; i++)
		for (j = 0; j < p; j++) {
		    tmp = 0.0;
		    for (k = 0; k < p; k++)
			tmp += T[i + p * k] * mm[k + p * j];
		    L[i + p * j] = tmp;
		}
	    for (i = 0; i < p; i++) {
		tmp = 0.0;
		for (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 (i = 0; i < p; i++)
		    for (j = 0; j < p; j++) {
			tmp = 0.0;
			for (k = 0; k < p; k++)
			    tmp += L[k + p * i] * Nt[l + 1 + n*k + n*p*j];
			mm[i + p * j] = tmp;
		    }
		for (i = 0; i < p; i++)
		    for (j = 0; j < p; j++) {
			tmp = 0.0;
			for (k = 0; k < p; k++)
			    tmp += mm[i + p * k] * L[k + p * j];
			Nt[l + n*i + n*p*j] += tmp;
		    }
	    }
	}

	for (i = 0; i < p; i++) {
	    tmp = 0.0;
	    for (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 (l = 0; l < n; l++) {
	    for (i = 0; i < p; i++)
		for (j = 0; j < p; j++) {
		    tmp = 0.0;
		    for (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 (i = 0; i < p; i++)
		for (j = 0; j < p; j++) {
		    tmp = Pt[l + n*i + n*p*j];
		    for (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 sZ, SEXP sa0, SEXP sP0, SEXP sT, SEXP sV,
	   SEXP sh, SEXP fast)
{
    SEXP res, forecasts, se;
    int  n = asReal(nahead), p = LENGTH(sa0);
    double *Z = REAL(sZ), *a = REAL(sa0), *P = REAL(sP0), *T = REAL(sT),
	*V = REAL(sV), h = asReal(sh);
    int i, j, k, l;
    double fc, tmp, *mm, *anew, *Pnew;

    /* It would be better to check types before using LENGTH and REAL
       on these, but should still work this way.  LT */
    if (TYPEOF(sZ) != REALSXP ||
	TYPEOF(sa0) != REALSXP || TYPEOF(sP0) != REALSXP ||
	TYPEOF(sT) != REALSXP || TYPEOF(sV) != REALSXP)
	error(_("invalid argument type"));

    anew = (double *) R_alloc(p, sizeof(double));
    Pnew = (double *) R_alloc(p * p, sizeof(double));
    mm = (double *) R_alloc(p * p, sizeof(double));
    PROTECT(res = allocVector(VECSXP, 2));
    SET_VECTOR_ELT(res, 0, forecasts = allocVector(REALSXP, n));
    SET_VECTOR_ELT(res, 1, se = allocVector(REALSXP, n));
    if (!LOGICAL(fast)[0]){
	PROTECT(sa0=duplicate(sa0));
	a=REAL(sa0);
	PROTECT(sP0=duplicate(sP0));
	P=REAL(sP0);
    }
    for (l = 0; l < n; l++) {
	fc = 0.0;
	for (i = 0; i < p; i++) {
	    tmp = 0.0;
	    for (k = 0; k < p; k++)
		tmp += T[i + p * k] * a[k];
	    anew[i] = tmp;
	    fc += tmp * Z[i];
	}
	for (i = 0; i < p; i++)
	    a[i] = anew[i];
	REAL(forecasts)[l] = fc;

	for (i = 0; i < p; i++)
	    for (j = 0; j < p; j++) {
		tmp = 0.0;
		for (k = 0; k < p; k++)
		    tmp += T[i + p * k] * P[k + p * j];
		mm[i + p * j] = tmp;
	    }
	for (i = 0; i < p; i++)
	    for (j = 0; j < p; j++) {
		tmp = V[i + p * j];
		for (k = 0; k < p; k++)
		    tmp += mm[i + p * k] * T[j + p * k];
		Pnew[i + p * j] = tmp;
	    }
	tmp = h;
	for (i = 0; i < p; i++)
	    for (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;
    }
    UNPROTECT(1);
    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],
	i, v, n = LENGTH(sin);
    double *params, *in = REAL(sin);
    SEXP res = allocVector(REALSXP, n);

    params = REAL(res);
    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);
    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;
}


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],
	i, j, v, n = LENGTH(in);
    SEXP y = allocMatrix(REALSXP, n, n);
    double *raw = REAL(in), *A = REAL(y), w1[100], w2[100], w3[100];

    for(i = 0; i < n; i++)
	for(j = 0; j < n; j++)
	    A[i + j*n] = (i == j);
    if(mp > 0) {
	for(i = 0; i < mp; i++) w1[i] = raw[i];
	partrans(mp, w1, w2);
	for(i = 0; i < mp; i++) {
	    w1[i] += eps;
	    partrans(mp, w1, w3);
	    for(j = 0; j < mp; j++) A[i + j*n] = (w3[j] - w2[j])/eps;
	    w1[i] -= eps;
	}
    }
    if(msp > 0) {
	v = mp + mq;
	for(i = 0; i < msp; i++) w1[i] = raw[i + v];
	partrans(msp, w1, w2);
	for(i = 0; i < msp; i++) {
	    w1[i] += eps;
	    partrans(msp, w1, w3);
	    for(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 sPhi, SEXP sTheta, SEXP sDelta,
	   SEXP sa, SEXP sP, SEXP sPn, SEXP sUP, SEXP giveResid)
{
    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, resid, gain, tmp, vi, *anew, *mm = NULL,
	*M;
    int i, j, k, l, 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 (l = 0; l < n; l++) {
	for (i = 0; i < r; i++) {
	    tmp = (i < r - 1) ? a[i + 1] : 0.0;
	    if (i < p) tmp += phi[i] * a[0];
	    anew[i] = tmp;
	}
	if (d > 0) {
	    for (i = r + 1; i < rd; i++) anew[i] = a[i - 1];
	    tmp = a[0];
	    for (i = 0; i < d; i++) tmp += delta[i] * a[r + i];
	    anew[r] = tmp;
	}
	if (l > asInteger(sUP)) {
	    if (d == 0) {
		for (i = 0; i < r; i++) {
		    vi = 0.0;
		    if (i == 0) vi = 1.0; else if (i - 1 < q) vi = theta[i - 1];
		    for (j = 0; j < r; j++) {
			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 (i = 0; i < r; i++)
		    for (j = 0; j < rd; j++) {
			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 (j = 0; j < rd; j++) {
		    tmp = P[rd * j];
		    for (k = 0; k < d; k++)
			tmp += delta[k] * P[r + k + rd * j];
		    mm[r + rd * j] = tmp;
		}
		for (i = 1; i < d; i++)
		    for (j = 0; j < rd; j++)
			mm[r + i + rd * j] = P[r + i - 1 + rd * j];

		/* Pnew = mmT' */
		for (i = 0; i < r; i++)
		    for (j = 0; j < rd; j++) {
			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 (j = 0; j < rd; j++) {
		    tmp = mm[j];
		    for (k = 0; k < d; k++)
			tmp += delta[k] * mm[rd * (r + k) + j];
		    Pnew[rd * r + j] = tmp;
		}
		for (i = 1; i < d; i++)
		    for (j = 0; j < rd; j++)
			Pnew[rd * (r + i) + j] = mm[rd * (r + i - 1) + j];
		/* Pnew <- Pnew + (1 theta) %o% (1 theta) */
		for (i = 0; i <= q; i++) {
		    vi = (i == 0) ? 1. : theta[i - 1];
		    for (j = 0; j <= q; j++)
			Pnew[i + rd * j] += vi * ((j == 0) ? 1. : theta[j - 1]);
		}
	    }
	}
	if (!ISNAN(y[l])) {
	    resid = y[l] - anew[0];
	    for (i = 0; i < d; i++)
		resid -= delta[i] * anew[r + i];

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

	    gain = M[0];
	    for (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 (i = 0; i < rd; i++)
		a[i] = anew[i] + M[i] * resid / gain;
	    for (i = 0; i < rd; i++)
		for (j = 0; j < rd; j++)
		    P[i + j * rd] = Pnew[i + j * rd] - M[i] * M[j] / gain;
	} else {
	    for (i = 0; i < rd; i++) a[i] = anew[i];
	    for (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 l, i, j, ns, nu = 0;
    Rboolean useResid = asLogical(giveResid);

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

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

    for (l = ncond; l < n; l++) {
	tmp = w[l];
	for (j = 0; j < p; j++) tmp -= phi[j] * w[l - j - 1];
	for (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 i, j, 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 (i = 0; i < nab; i++) rab[i] = 0.0;
    for (i = 0; i < na; i++)
	for (j = 0; j < nb; j++)
	    rab[i + j] += ra[i] * rb[j];
    UNPROTECT(3);
    return (ab);
}


inline int triangular(int row, int col, int N)
{
   if (row <= col)
      return row*N - (row-1)*((row-1) + 1)/2 + col - row;
   else if (col<row)
      return col*N - (col-1)*((col-1) + 1)/2 + row - col;
}

inline void dot(double *A, int rowsA, int colsA, double *B, int colsB, double *C) {
    // C <- dot(A, B)
    int i, j, k, ind;
    for (i = 0; i < rowsA; i++) {
        for (j = 0; j < colsB; j++) {
	    ind = i * colsB + j;
	    C[ind] = 0.0;
	    for (k = 0; k < colsA; k++) {
	        C[ind] += A[i * colsA + k] * B[k * colsB + j];
	    }
	}
    }
}

inline void dotT(double *A, int rowsA, int colsA, double *B, int colsB, double *C) {
    // C <- dot(A, B.T)
    int i, j, k, ind;
    for (i = 0; i < rowsA; i++) {
        for (j = 0; j < colsB; j++) {
	    ind = i * colsB + j;
	    C[ind] = 0.0;
	    for (k = 0; k < colsA; k++) {
	        C[ind] += A[i * colsA + k] * B[j * colsA + k];
	    }
	}
    }
}

inline void Tdot(double *T, double *B, int r, double *C) {
    // C <- dot(T, B)
    int i, j, k, ind;
    for (i = 0; i < r; i++) {
        for (j = 0; j < r; j++) {
	    ind = i * r + j;
	    k = 0;
	    C[ind] = T[i * r + k] * B[k * r + j];
	    k = i + 1;
	    if (k < r) {
		C[ind] += T[i * r + k] * B[k * r + j];
	    }
	}
    }
}

inline void TdotTplusV(double *A, double *T, int r, double *C, double *D) {
    // C <- dot(A, T') + V
    int i, j, k, ind;
    for (i = 0; i < r; i++) {
        for (j = 0; j < r; j++) {
	    ind = i * r + j;
	    k = 0;
	    C[ind] = D[ind] + A[i * r + k] * T[j * r + k];
	    k = j + 1;
	    if (k < r) {
		C[ind] += A[i * r + k] * T[j * r + k];
	    }
	}
    }
}

inline void add(double *A, int rowsA, int colsA, double *B) {
    // A <- add(A, B)
    int i, j, ind;
    for (i = 0; i < rowsA; i++) {
        for (j = 0; j < colsA; j++) {
	    ind = i * colsA + j;
	    A[ind] += B[ind];
	}
    }
}

inline void eye(double *A, int rowsA, int colsA) {
    int i, j, k, ind;
    for (i = 0; i < rowsA; i++) {
        for (j = 0; j < colsA; j++) {
	    if (i == j) A[i * rowsA + j] = 1.0;
	    else A[i * rowsA + j] = 0.0;
	}
    }
}

/*
 * Based on code from subroutine STARMA in AS154
 * http://lib.stat.cmu.edu/apstat/154
 * */
SEXP getQ0(SEXP sPhi, SEXP sTheta) //, SEXP sToll)
{
    SEXP res;
    int  p = LENGTH(sPhi), q = LENGTH(sTheta);
    double *phi = REAL(sPhi), *theta = REAL(sTheta);
    //double toll = asReal(sToll);

    double *P, *V, *R, *T, *Pnew;
    int r = max(p, q + 1);
    size_t np = r * (r + 1) / 2;
    int indi, indj, indn;
    double phii, phij, ynext, bi, vi, vj, tk, tl;
    int   i, j, k, l, ithisr, ind, npr, ind1, ind2, npr1, im, jm;

    R = (double *) R_alloc(r, sizeof(double));
    /*
     *     /     1      \
     *     |  theta[0]  |
     * R = |  theta[1]  |
     *     |     .      |
     *     \ theta[r-1] /
     *
     * */
    R[0] = 1.0;
    for (i = 1; i < r; i++) {
        if (i <= q) R[i] = theta[i - 1];
	else R[i] = 0.0;
    }

    T = (double *) R_alloc(r*r, sizeof(double));
    /*
     *     / phi[0] 1 0 . 0 \
     *     | phi[1] 0 1 0 0 |
     * T = |   .        .   |
     *     |   .          1 |
     *     \ phi[r] 0 . . 0 /
     *
     * */
    for (i = 0; i < r; i++) {
        for (j = 0; j < r; j++) {
	    if (j == 0 && i < p) T[i * r + j] = phi[i];
	    else if ((i+1) == j) T[i * r + j] = 1.0;
	    else T[i * r + j] = 0.0;
	}
    }

    V = (double *) R_alloc(r*r, sizeof(double));
    dot(R, r, 1, R, r, V);

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

    if (r == 1) { /* p == 1, q == 0 */
	P[0] = 1.0 / (1.0 - phi[0] * phi[0]);
	UNPROTECT(1);
	return res;
    }
    /* if (p > 0) { */
    /* Iteratively solve: P0 = T*P0*T' + R*R' */
	eye(P, r, r);
	for (ind = 0; ind < (2 * r); ind++) {
	    Tdot(T, P, r, Pnew);
	    TdotTplusV(Pnew, T, r, P, V);
	}
    /* } 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; // res is the matrix view of P
}

/* vim: set ts=8: */

Hi, is there any reason why parameter transformation is limited to just 100 parameters?
I was trying to fit model to hourly data with weekly seasonality via arfimafit function in rugarch package.
My model is of order 24 * 7 = 168 with most parameters fixed to 0.

I hardly remember anything about this code, but that 100 seems to me just a constant number to avoid dynamic allocation, set it to 200 and it should go @danaszuokas

Do you remember what references you used for this specific script?

@git314 beside what I left as comment (https://gist.github.com/naufraghi/911292/revisions shows the modifications I did on the sources) I forgot most details, I think the "iterative solution" here https://gist.github.com/naufraghi/911292#file-arima-c-L935 is not verified, I simply noted that in most cases that worked and avoided a huge allocation here https://gist.github.com/naufraghi/911292/revisions#diff-e2dc2642f747203273aeaf8d4eb2dc097f0c045dc863aad5cbbd0a67f7ee6b6eL842. The paper that used this code should be this one: https://ieeexplore.ieee.org/document/6482260