/* nag_tsa_transf_filter (g13bbc) 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, ii, ij, iqxd, j, k, n, nb, ni, npar, nparx, nx;
Integer nser, npara, tdxxy, tdmrx, ldparx, tdparx;
Integer exit_status = 0, idd = 0, ny = 0;
/* Arrays */
double *b = 0, *fsd = 0, *fva = 0, *par = 0, *parx = 0;
double *x = 0, *y = 0, *rms = 0, *parxx = 0;
Integer mr[10], mrx[7], *mrxx = 0;
Nag_TransfOrder transfj, transfv;
Nag_ArimaOrder arimaj, arimas;
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_transf_filter (g13bbc) Example Program Results\n");
/* Skip heading in data file */
scanf("%*[^\n] ");
scanf("%ld%*[^\n] ", &nx);
printf("\n");
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 <= 3; ++i)
scanf("%ld", &mr[i-1]);
scanf("%*[^\n] ");
transfv.nag_b = mr[0];
transfv.nag_q = mr[1];
transfv.nag_p = mr[2];
npar = mr[1] + mr[2] + 1;
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];
if (iqxd != 0)
{
if (!(fsd = NAG_ALLOC(iqxd, double)) ||
!(fva = NAG_ALLOC(iqxd, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
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;
}
}
/* Calculate series length */
ny = nx + iqxd;
/* Allocate array y */
if (!(y = NAG_ALLOC(ny, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Move backforecasts to start of y array */
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 > 215)
goto END;
y[j-1] = x[k-1];
++j;
--k;
}
}
/* Move ARIMA for series into mr */
for (i = 1; i <= 7; ++i)
mr[i+2] = mrx[i-1];
arimas.p = mr[3];
arimas.d = mr[4];
arimas.q = mr[5];
arimas.bigp = mr[6];
arimas.bigd = mr[7];
arimas.bigq = mr[8];
arimas.s = mr[9];
/* 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_transf_filter (g13bbc) */
nb = ny + MAX(mr[0] + mr[1], mr[2]);
/* 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_transf_filter (g13bbc) */
/* nag_tsa_transf_filter (g13bbc).
* Multivariate time series, filtering by a transfer
* function model
*/
nag_tsa_transf_filter(y, ny, &transfv, &arimas, par, npar, cy, b, nb,
&fail);
if (fail.code != NE_NOERROR)
{
printf(
"Error from nag_tsa_transf_filter (g13bbc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
printf(" Original Filtered\n");
printf(" Backforecasts y-series series\n");
if (iqxd != 0)
{
ij = -iqxd;
for (i = 1; i <= iqxd; ++i)
{
printf("%8ld%17.1f%16.1f\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("%10.1f", 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;
}