1 /* R : A Computer Language for Statistical Data Analysis
2 *
3 * Copyright (C) 1999-2016 The R Core Team
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, a copy is available at
17 * https://www.R-project.org/Licenses/.
18 */
19
20 #ifdef HAVE_CONFIG_H
21 # include <config.h>
22 #endif
23
24 #include <R.h>
25 #include "ts.h"
26
27 #ifndef max
28 #define max(a,b) ((a < b)?(b):(a))
29 #endif
30 #ifndef min
31 #define min(a,b) ((a > b)?(b):(a))
32 #endif
33
34
35 /* Internal */
36 static void partrans(int np, double *raw, double *new);
37 static void dotrans(Starma G, double *raw, double *new, int trans);
38
39
40 /* cor is the autocorrelations starting from 0 lag*/
uni_pacf(double * cor,double * p,int nlag)41 static void uni_pacf(double *cor, double *p, int nlag)
42 {
43 double a, b, c, *v, *w;
44
45 v = (double*) R_alloc(nlag, sizeof(double));
46 w = (double*) R_alloc(nlag, sizeof(double));
47 w[0] = p[0] = cor[1];
48 for(int ll = 1; ll < nlag; ll++) {
49 a = cor[ll+1];
50 b = 1.0;
51 for(int i = 0; i < ll; i++) {
52 a -= w[i] * cor[ll - i];
53 b -= w[i] * cor[i + 1];
54 }
55 p[ll] = c = a/b;
56 if(ll+1 == nlag) break;
57 w[ll] = c;
58 for(int i = 0; i < ll; i++)
59 v[ll-i-1] = w[i];
60 for(int i = 0; i < ll; i++)
61 w[i] -= c*v[i];
62 }
63 }
64
pacf1(SEXP acf,SEXP lmax)65 SEXP pacf1(SEXP acf, SEXP lmax)
66 {
67 int lagmax = asInteger(lmax);
68 acf = PROTECT(coerceVector(acf, REALSXP));
69 SEXP ans = PROTECT(allocVector(REALSXP, lagmax));
70 uni_pacf(REAL(acf), REAL(ans), lagmax);
71 SEXP d = PROTECT(allocVector(INTSXP, 3));
72 INTEGER(d)[0] = lagmax;
73 INTEGER(d)[1] = INTEGER(d)[2] = 1;
74 setAttrib(ans, R_DimSymbol, d);
75 UNPROTECT(3);
76 return ans;
77 }
78
79
80 /* Use an external reference to store the structure we keep allocated
81 memory in */
82 static SEXP Starma_tag;
83
84 #define GET_STARMA \
85 Starma G; \
86 if (TYPEOF(pG) != EXTPTRSXP || R_ExternalPtrTag(pG) != Starma_tag) \
87 error(_("bad Starma struct"));\
88 G = (Starma) R_ExternalPtrAddr(pG)
89
setup_starma(SEXP na,SEXP x,SEXP pn,SEXP xreg,SEXP pm,SEXP dt,SEXP ptrans,SEXP sncond)90 SEXP setup_starma(SEXP na, SEXP x, SEXP pn, SEXP xreg, SEXP pm,
91 SEXP dt, SEXP ptrans, SEXP sncond)
92 {
93 Starma G;
94 int i, n, m, ip, iq, ir, np;
95 SEXP res;
96 double *rx = REAL(x), *rxreg = REAL(xreg);
97
98 G = Calloc(1, starma_struct);
99 G->mp = INTEGER(na)[0];
100 G->mq = INTEGER(na)[1];
101 G->msp = INTEGER(na)[2];
102 G->msq = INTEGER(na)[3];
103 G->ns = INTEGER(na)[4];
104 G->n = n = asInteger(pn);
105 G->ncond = asInteger(sncond);
106 G->m = m = asInteger(pm);
107 G->params = Calloc(G->mp + G->mq + G->msp + G->msq + G->m, double);
108 G->p = ip = G->ns*G->msp + G->mp;
109 G->q = iq = G->ns*G->msq + G->mq;
110 G->r = ir = max(ip, iq + 1);
111 G->np = np = (ir*(ir + 1))/2;
112 G->nrbar = max(1, np*(np - 1)/2);
113 G->trans = asInteger(ptrans);
114 G->a = Calloc(ir, double);
115 G->P = Calloc(np, double);
116 G->V = Calloc(np, double);
117 G->thetab = Calloc(np, double);
118 G->xnext = Calloc(np, double);
119 G->xrow = Calloc(np, double);
120 G->rbar = Calloc(G->nrbar, double);
121 G->w = Calloc(n, double);
122 G->wkeep = Calloc(n, double);
123 G->resid = Calloc(n, double);
124 G->phi = Calloc(ir, double);
125 G->theta = Calloc(ir, double);
126 G->reg = Calloc(1 + n*m, double); /* AIX can't calloc 0 items */
127 G->delta = asReal(dt);
128 for(i = 0; i < n; i++) G->w[i] = G->wkeep[i] = rx[i];
129 for(i = 0; i < n*m; i++) G->reg[i] = rxreg[i];
130 Starma_tag = install("STARMA_TAG");
131 res = R_MakeExternalPtr(G, Starma_tag, R_NilValue);
132 return res;
133 }
134
free_starma(SEXP pG)135 SEXP free_starma(SEXP pG)
136 {
137 GET_STARMA;
138
139 Free(G->params); Free(G->a); Free(G->P); Free(G->V); Free(G->thetab);
140 Free(G->xnext); Free(G->xrow); Free(G->rbar);
141 Free(G->w); Free(G->wkeep); Free(G->resid); Free(G->phi); Free(G->theta);
142 Free(G->reg); Free(G);
143 return R_NilValue;
144 }
145
Starma_method(SEXP pG,SEXP method)146 SEXP Starma_method(SEXP pG, SEXP method)
147 {
148 GET_STARMA;
149
150 G->method = asInteger(method);
151 return R_NilValue;
152 }
153
Dotrans(SEXP pG,SEXP x)154 SEXP Dotrans(SEXP pG, SEXP x)
155 {
156 SEXP y = allocVector(REALSXP, LENGTH(x));
157 GET_STARMA;
158
159 dotrans(G, REAL(x), REAL(y), 1);
160 return y;
161 }
162
set_trans(SEXP pG,SEXP ptrans)163 SEXP set_trans(SEXP pG, SEXP ptrans)
164 {
165 GET_STARMA;
166
167 G->trans = asInteger(ptrans);
168 return R_NilValue;
169 }
170
arma0fa(SEXP pG,SEXP inparams)171 SEXP arma0fa(SEXP pG, SEXP inparams)
172 {
173 int i, j, ifault = 0, it, streg;
174 double sumlog, ssq, tmp, ans;
175
176 GET_STARMA;
177 dotrans(G, REAL(inparams), G->params, G->trans);
178
179 if(G->ns > 0) {
180 /* expand out seasonal ARMA models */
181 for(i = 0; i < G->mp; i++) G->phi[i] = G->params[i];
182 for(i = 0; i < G->mq; i++) G->theta[i] = G->params[i + G->mp];
183 for(i = G->mp; i < G->p; i++) G->phi[i] = 0.0;
184 for(i = G->mq; i < G->q; i++) G->theta[i] = 0.0;
185 for(j = 0; j < G->msp; j++) {
186 G->phi[(j + 1)*G->ns - 1] += G->params[j + G->mp + G->mq];
187 for(i = 0; i < G->mp; i++)
188 G->phi[(j + 1)*G->ns + i] -= G->params[i]*
189 G->params[j + G->mp + G->mq];
190 }
191 for(j = 0; j < G->msq; j++) {
192 G->theta[(j + 1)*G->ns - 1] +=
193 G->params[j + G->mp + G->mq + G->msp];
194 for(i = 0; i < G->mq; i++)
195 G->theta[(j + 1)*G->ns + i] += G->params[i + G->mp]*
196 G->params[j + G->mp + G->mq + G->msp];
197 }
198 } else {
199 for(i = 0; i < G->mp; i++) G->phi[i] = G->params[i];
200 for(i = 0; i < G->mq; i++) G->theta[i] = G->params[i + G->mp];
201 }
202
203 streg = G->mp + G->mq + G->msp + G->msq;
204 if(G->m > 0) {
205 for(i = 0; i < G->n; i++) {
206 tmp = G->wkeep[i];
207 for(j = 0; j < G->m; j++)
208 tmp -= G->reg[i + G->n*j] * G->params[streg + j];
209 G->w[i] = tmp;
210 }
211 }
212
213 if(G->method == 1) {
214 int p = G->mp + G->ns * G->msp, q = G->mq + G->ns * G->msq, nu = 0;
215 ssq = 0.0;
216 for(i = 0; i < G->ncond; i++) G->resid[i] = 0.0;
217 for(i = G->ncond; i < G->n; i++) {
218 tmp = G->w[i];
219 for(j = 0; j < min(i - G->ncond, p); j++)
220 tmp -= G->phi[j] * G->w[i - j - 1];
221 for(j = 0; j < min(i - G->ncond, q); j++)
222 tmp -= G->theta[j] * G->resid[i - j - 1];
223 G->resid[i] = tmp;
224 if(!ISNAN(tmp)) {
225 nu++;
226 ssq += tmp * tmp;
227 }
228 }
229 G->s2 = ssq/(double)(nu);
230 ans = 0.5 * log(G->s2);
231 } else {
232 starma(G, &ifault);
233 if(ifault) error(_("starma error code %d"), ifault);
234 sumlog = 0.0;
235 ssq = 0.0;
236 it = 0;
237 karma(G, &sumlog, &ssq, 1, &it);
238 G->s2 = ssq/(double)G->nused;
239 ans = 0.5*(log(ssq/(double)G->nused) + sumlog/(double)G->nused);
240 }
241 return ScalarReal(ans);
242 }
243
get_s2(SEXP pG)244 SEXP get_s2(SEXP pG)
245 {
246 GET_STARMA;
247 return ScalarReal(G->s2);
248 }
249
get_resid(SEXP pG)250 SEXP get_resid(SEXP pG)
251 {
252 SEXP res;
253 int i;
254 double *rres;
255 GET_STARMA;
256
257 res = allocVector(REALSXP, G->n);
258 rres = REAL(res);
259 for(i = 0; i < G->n; i++) rres[i] = G->resid[i];
260 return res;
261 }
262
arma0_kfore(SEXP pG,SEXP pd,SEXP psd,SEXP nahead)263 SEXP arma0_kfore(SEXP pG, SEXP pd, SEXP psd, SEXP nahead)
264 {
265 int dd = asInteger(pd);
266 int d, il = asInteger(nahead), ifault = 0, i, j;
267 double *del, *del2;
268 SEXP res, x, var;
269 GET_STARMA;
270
271 PROTECT(res = allocVector(VECSXP, 2));
272 SET_VECTOR_ELT(res, 0, x = allocVector(REALSXP, il));
273 SET_VECTOR_ELT(res, 1, var = allocVector(REALSXP, il));
274
275 d = dd + G->ns * asInteger(psd);
276
277 del = (double *) R_alloc(d + 1, sizeof(double));
278 del2 = (double *) R_alloc(d + 1, sizeof(double));
279 del[0] = 1;
280 for(i = 1; i <= d; i++) del[i] = 0;
281 for (j = 0; j < dd; j++) {
282 for(i = 0; i <= d; i++) del2[i] = del[i];
283 for(i = 0; i <= d - 1; i++) del[i+1] -= del2[i];
284 }
285 for (j = 0; j < asInteger(psd); j++) {
286 for(i = 0; i <= d; i++) del2[i] = del[i];
287 for(i = 0; i <= d - G->ns; i++) del[i + G->ns] -= del2[i];
288 }
289 for(i = 1; i <= d; i++) del[i] *= -1;
290
291
292 forkal(G, d, il, del + 1, REAL(x), REAL(var), &ifault);
293 if(ifault) error(_("forkal error code %d"), ifault);
294 UNPROTECT(1);
295 return res;
296 }
297
artoma(int p,double * phi,double * psi,int npsi)298 static void artoma(int p, double *phi, double *psi, int npsi)
299 {
300 int i, j;
301
302 for(i = 0; i < p; i++) psi[i] = phi[i];
303 for(i = p; i < npsi; i++) psi[i] = 0.0;
304 for(i = 0; i < npsi - p - 1; i++)
305 for(j = 0; j < p; j++) psi[i + j + 1] += phi[j]*psi[i];
306 }
307
ar2ma(SEXP ar,SEXP npsi)308 SEXP ar2ma(SEXP ar, SEXP npsi)
309 {
310 ar = PROTECT(coerceVector(ar, REALSXP));
311 int p = LENGTH(ar), ns = asInteger(npsi), ns1 = ns + p + 1;
312 SEXP psi = PROTECT(allocVector(REALSXP, ns1));
313 artoma(p, REAL(ar), REAL(psi), ns1);
314 SEXP ans = lengthgets(psi, ns);
315 UNPROTECT(2);
316 return ans;
317 }
318
partrans(int p,double * raw,double * new)319 static void partrans(int p, double *raw, double *new)
320 {
321 int j, k;
322 double a, work[100];
323
324 if(p > 100) error(_("can only transform 100 pars in arima0"));
325
326 /* Step one: map (-Inf, Inf) to (-1, 1) via tanh
327 The parameters are now the pacf phi_{kk} */
328 for(j = 0; j < p; j++) work[j] = new[j] = tanh(raw[j]);
329 /* Step two: run the Durbin-Levinson recursions to find phi_{j.},
330 j = 2, ..., p and phi_{p.} are the autoregression coefficients */
331 for(j = 1; j < p; j++) {
332 a = new[j];
333 for(k = 0; k < j; k++)
334 work[k] -= a * new[j - k - 1];
335 for(k = 0; k < j; k++) new[k] = work[k];
336 }
337 }
338
dotrans(Starma G,double * raw,double * new,int trans)339 static void dotrans(Starma G, double *raw, double *new, int trans)
340 {
341 int i, v, n = G->mp + G->mq + G->msp + G->msq + G->m;
342
343 for(i = 0; i < n; i++) new[i] = raw[i];
344 if(trans) {
345 partrans(G->mp, raw, new);
346 v = G->mp;
347 partrans(G->mq, raw + v, new + v);
348 v += G->mq;
349 partrans(G->msp, raw + v, new + v);
350 v += G->msp;
351 partrans(G->msq, raw + v, new + v);
352 }
353 }
354
355 #if !defined(atanh) && defined(HAVE_DECL_ATANH) && !HAVE_DECL_ATANH
356 extern double atanh(double x);
357 #endif
invpartrans(int p,double * phi,double * new)358 static void invpartrans(int p, double *phi, double *new)
359 {
360 int j, k;
361 double a, work[100];
362
363 if(p > 100) error(_("can only transform 100 pars in arima0"));
364
365 for(j = 0; j < p; j++) work[j] = new[j] = phi[j];
366 /* Run the Durbin-Levinson recursions backwards
367 to find the PACF phi_{j.} from the autoregression coefficients */
368 for(j = p - 1; j > 0; j--) {
369 a = new[j];
370 for(k = 0; k < j; k++)
371 work[k] = (new[k] + a * new[j - k - 1]) / (1 - a * a);
372 for(k = 0; k < j; k++) new[k] = work[k];
373 }
374 for(j = 0; j < p; j++) new[j] = atanh(new[j]);
375 }
376
Invtrans(SEXP pG,SEXP x)377 SEXP Invtrans(SEXP pG, SEXP x)
378 {
379 SEXP y = allocVector(REALSXP, LENGTH(x));
380 int i, v, n;
381 double *raw = REAL(x), *new = REAL(y);
382 GET_STARMA;
383
384 n = G->mp + G->mq + G->msp + G->msq;
385
386 v = 0;
387 invpartrans(G->mp, raw + v, new + v);
388 v += G->mp;
389 invpartrans(G->mq, raw + v, new + v);
390 v += G->mq;
391 invpartrans(G->msp, raw + v, new + v);
392 v += G->msp;
393 invpartrans(G->msq, raw + v, new + v);
394 for(i = n; i < n + G->m; i++) new[i] = raw[i];
395 return y;
396 }
397
398 #define eps 1e-3
Gradtrans(SEXP pG,SEXP x)399 SEXP Gradtrans(SEXP pG, SEXP x)
400 {
401 SEXP y = allocMatrix(REALSXP, LENGTH(x), LENGTH(x));
402 int i, j, v, n;
403 double *raw = REAL(x), *A = REAL(y), w1[100], w2[100], w3[100];
404 GET_STARMA;
405
406 n = G->mp + G->mq + G->msp + G->msq + G->m;
407 for(i = 0; i < n; i++)
408 for(j = 0; j < n; j++)
409 A[i + j*n] = (i == j);
410 if(G->mp > 0) {
411 for(i = 0; i < G->mp; i++) w1[i] = raw[i];
412 partrans(G->mp, w1, w2);
413 for(i = 0; i < G->mp; i++) {
414 w1[i] += eps;
415 partrans(G->mp, w1, w3);
416 for(j = 0; j < G->mp; j++) A[i + j*n] = (w3[j] - w2[j])/eps;
417 w1[i] -= eps;
418 }
419 }
420 if(G->mq > 0) {
421 v = G->mp;
422 for(i = 0; i < G->mq; i++) w1[i] = raw[i + v];
423 partrans(G->mq, w1, w2);
424 for(i = 0; i < G->mq; i++) {
425 w1[i] += eps;
426 partrans(G->mq, w1, w3);
427 for(j = 0; j < G->mq; j++) A[i + v + j*n] = (w3[j] - w2[j])/eps;
428 w1[i] -= eps;
429 }
430 }
431 if(G->msp > 0) {
432 v = G->mp + G->mq;
433 for(i = 0; i < G->msp; i++) w1[i] = raw[i + v];
434 partrans(G->msp, w1, w2);
435 for(i = 0; i < G->msp; i++) {
436 w1[i] += eps;
437 partrans(G->msp, w1, w3);
438 for(j = 0; j < G->msp; j++)
439 A[i + v + (j+v)*n] = (w3[j] - w2[j])/eps;
440 w1[i] -= eps;
441 }
442 }
443 if(G->msq > 0) {
444 v = G->mp + G->mq + G->msp;
445 for(i = 0; i < G->msq; i++) w1[i] = raw[i + v];
446 partrans(G->msq, w1, w2);
447 for(i = 0; i < G->msq; i++) {
448 w1[i] += eps;
449 partrans(G->msq, w1, w3);
450 for(j = 0; j < G->msq; j++)
451 A[i + v + (j+v)*n] = (w3[j] - w2[j])/eps;
452 w1[i] -= eps;
453 }
454 }
455 return y;
456 }
457
458 SEXP
ARMAtoMA(SEXP ar,SEXP ma,SEXP lag_max)459 ARMAtoMA(SEXP ar, SEXP ma, SEXP lag_max)
460 {
461 int i, j, p = LENGTH(ar), q = LENGTH(ma), m = asInteger(lag_max);
462 double *phi = REAL(ar), *theta = REAL(ma), *psi, tmp;
463 SEXP res;
464
465 if(m <= 0 || m == NA_INTEGER)
466 error(_("invalid value of lag.max"));
467 PROTECT(res = allocVector(REALSXP, m));
468 psi = REAL(res);
469 for(i = 0; i < m; i++) {
470 tmp = (i < q) ? theta[i] : 0.0;
471 for(j = 0; j < min(i+1, p); j++)
472 tmp += phi[j] * ((i-j-1 >= 0) ? psi[i-j-1] : 1.0);
473 psi[i] = tmp;
474 }
475 UNPROTECT(1);
476 return res;
477 }
478