/* nag_tsa_arma_filter (g13bac) Example Program.
*
* Copyright 2014 Numerical Algorithms Group.
*
* Mark 7, 2002.
* Mark 7b revised, 2004.
*/
#include <stdio.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagg13.h>
int main(void)
{
/* Scalars */
double a1, a2, cx, cy;
Integer i, idd, ii, ij, iqxd,
j, k, n, nb, ni, npar, nparx, nx, ny,
nser, npara, tdxxy, tdmrx, ldparx, tdparx;
Integer exit_status = 0;
/* Arrays */
double *b = 0, *fsd = 0, *fva = 0, *par = 0, *parx = 0,
*x = 0, *y = 0, *rms = 0, *parxx = 0;
Integer mr[14], mrx[7], *mrxx = 0;
Nag_ArimaOrder arimaj, arimaf, arimav;
Nag_TransfOrder transfj;
Nag_G13_Opt options;
NagError fail;
INIT_FAIL(fail);
exit_status = 0;
/* Initialise the options structure used by nag_tsa_multi_inp_model_forecast
* (g13bjc) */
/* nag_tsa_options_init (g13bxc).
* Initialization function for option setting
*/
nag_tsa_options_init(&options);
printf("nag_tsa_arma_filter (g13bac) Example Program Results\n");
/* Skip heading in data file */
scanf("%*[^\n] ");
scanf("%ld%*[^\n] ", &nx);
if (nx > 0)
{
/* Allocate array x */
if (!(x = NAG_ALLOC(nx+2, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
for (i = 1; i <= nx; ++i)
scanf("%lf", &x[i-1]);
scanf("%*[^\n] ");
/* Read univariate Arima for series */
for (i = 1; i <= 7; ++i)
scanf("%ld", &mrx[i-1]);
scanf("%*[^\n] ");
scanf("%lf%*[^\n] ", &cx);
nparx = mrx[0] + mrx[2] + mrx[3] + mrx[5];
arimaj.p = mrx[0];
arimaj.d = mrx[1];
arimaj.q = mrx[2];
arimaj.bigp = mrx[3];
arimaj.bigd = mrx[4];
arimaj.bigq = mrx[5];
arimaj.s = mrx[6];
nser = 1;
if (nparx > 0)
{
/* Allocate array parx */
if (!(parx = NAG_ALLOC(nparx+nser, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
for (i = 1; i <= nparx; ++i)
scanf("%lf", &parx[i-1]);
scanf("%*[^\n] ");
/* Read model by which to filter series */
for (i = 1; i <= 7; ++i)
scanf("%ld", &mr[i-1]);
scanf("%*[^\n] ");
arimaf.p = mr[0];
arimaf.d = mr[1];
arimaf.q = mr[2];
arimaf.bigp = mr[3];
arimaf.bigd = mr[4];
arimaf.bigq = mr[5];
arimaf.s = mr[6];
npar = mr[0] + mr[2] + mr[3] + mr[5];
if (npar > 0)
{
/* Allocate array par */
if (!(par = NAG_ALLOC(npar + nparx, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
for (i = 1; i <= npar; ++i)
scanf("%lf", &par[i-1]);
scanf("%*[^\n] ");
/* Initially backforecast QY values */
/* (1) Reverse series in situ */
n = nx / 2;
ni = nx;
for (i = 1; i <= n; ++i)
{
a1 = x[i-1];
a2 = x[ni-1];
x[i-1] = a2;
x[ni-1] = a1;
--ni;
}
idd = mrx[1] + mrx[4];
/* (2) Possible sign reversal for ARIMA constant */
if (idd % 2 != 0)
cx = -cx;
/* (3) Calculate number of backforecasts required */
iqxd = mrx[2] + mrx[5] * mrx[6];
/* Calculate series length */
ny = nx + iqxd;
/* Allocate array y */
if (!(y = NAG_ALLOC(ny, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
if (iqxd != 0)
{
/* Allocate arrays fsd, fva and st. */
if (!(fsd = NAG_ALLOC(iqxd, double)) ||
!(fva = NAG_ALLOC(iqxd, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* (4) Set up parameter list for call to forecast
* routine g13bjc
*/
npara = nparx+nser;
parx[npara-1] = cx;
tdxxy = nser;
tdmrx = nser-1;
ldparx = nser-1;
tdparx = nser-1;
if (!(rms = NAG_ALLOC(nser, double)) ||
!(parxx = NAG_ALLOC(nser, double)) ||
!(mrxx = NAG_ALLOC(7*nser, Integer)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* nag_tsa_transf_orders (g13byc).
* Allocates memory to transfer function model orders
*/
nag_tsa_transf_orders(nser, &transfj, &fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_tsa_transf_orders (g13byc)"
".\n%s\n", fail.message);
exit_status = 1;
goto END;
}
rms[0] = 0;
transfj.nag_b = 0;
transfj.nag_q = 0;
transfj.nag_p = 0;
transfj.nag_r = 1;
for (i = 1; i <= 7; ++i)
mrxx[i-1] = 0;
parxx[0] = 0;
/* Tell nag_tsa_multi_inp_model_forecast (g13bjc) not to
* print parameters on entry */
options.list = Nag_FALSE;
/* nag_tsa_multi_inp_model_forecast (g13bjc).
* Forecasting function
*/
nag_tsa_multi_inp_model_forecast(&arimaj, nser, &transfj,
parx, npara, nx, iqxd, x,
tdxxy, rms, mrxx, tdmrx,
parxx, ldparx, tdparx,
fva, fsd, &options, &fail);
if (fail.code != NE_NOERROR)
{
printf(
"Error from nag_tsa_multi_inp_model_forecast "
"(g13bjc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
j = iqxd;
for (i = 1; i <= iqxd; ++i)
{
y[i-1] = fva[j-1];
--j;
}
/* Move series into y */
j = iqxd + 1;
k = nx;
for (i = 1; i <= nx; ++i)
{
if (j > 305)
goto END;
y[j-1] = x[k-1];
++j;
--k;
}
}
/* Move ARIMA for series into mr */
for (i = 1; i <= 7; ++i)
mr[i+6] = mrx[i-1];
arimav.p = mr[7];
arimav.d = mr[8];
arimav.q = mr[9];
arimav.bigp = mr[10];
arimav.bigd = mr[11];
arimav.bigq = mr[12];
arimav.s = mr[13];
/* Move parameters of ARIMA for y into par */
for (i = 1; i <= nparx; ++i)
par[npar+i-1] = parx[i-1];
npar += nparx;
/* Move constant and reset sign reversal */
cy = cx;
if (idd % 2 != 0)
cy = -cy;
/* Set parameters for call to filter routine
* nag_tsa_arma_filter (g13bac) */
nb = ny + MAX(mr[2] + mr[5] * mr[6],
mr[0] + mr[1] + (mr[3] + mr[4]) * mr[6]);
/* Allocate array b */
if (!(b = NAG_ALLOC(nb, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Filter series by call to nag_tsa_arma_filter (g13bac) */
/* nag_tsa_arma_filter (g13bac).
* Multivariate time series, filtering (pre-whitening) by an
* ARIMA model
*/
nag_tsa_arma_filter(y, ny, &arimaf, &arimav, par, npar, cy, b,
nb, &fail);
if (fail.code != NE_NOERROR)
{
printf(
"Error from nag_tsa_arma_filter (g13bac).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
printf("\n");
printf(" Original Filtered\n");
printf("Backforecasts y-series series\n");
if (iqxd != 0)
{
ij = -iqxd;
for (i = 1; i <= iqxd; ++i)
{
printf("%8ld%17.4f%15.4f\n", ij, y[i-1],
b[i-1]);
++ij;
}
}
printf("\n");
printf(" Filtered Filtered "
"Filtered Filtered\n");
printf(" series series "
" series series\n");
for (i = iqxd + 1; i <= ny; i += 4)
{
for (ii = i; ii <= MIN(ny, i+3); ++ii)
{
printf("%5ld", ii-iqxd);
printf("%9.4f ", b[ii-1]);
}
printf("\n");
}
}
}
}
END:
/* Free the options structure used by nag_tsa_multi_inp_model_forecast
* (g13bjc) */
/* nag_tsa_free (g13xzc).
* Freeing function for use with g13 option setting
*/
nag_tsa_free(&options);
NAG_FREE(b);
NAG_FREE(fsd);
NAG_FREE(fva);
NAG_FREE(par);
NAG_FREE(parx);
NAG_FREE(x);
NAG_FREE(y);
NAG_FREE(rms);
NAG_FREE(parxx);
NAG_FREE(mrxx);
return exit_status;
}