/* nag_tsa_transf_filter (g13bbc) Example Program. * * Copyright 2002 Numerical Algorithms Group. * * Mark 7, 2002. */ #include #include #include #include int main(void) { /* Scalars */ double a1, a2, cx, cy; Integer i, idd, ii, ij, iqxd, j, k, n, nb, ni, nmr, 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[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 g13bjc */ g13bxc(&options); Vprintf("g13bbc Example Program Results\n"); /* Skip heading in data file */ Vscanf("%*[^\n] "); Vscanf("%ld%*[^\n] ", &nx); Vprintf("\n"); if (nx > 0) { /* Allocate array x */ if ( !(x = NAG_ALLOC(nx+2, double))) { Vprintf("Allocation failure\n"); exit_status = -1; goto END; } for (i = 1; i <= nx; ++i) Vscanf("%lf", &x[i-1]); Vscanf("%*[^\n] "); /* Read univariate ARIMA for series */ for (i = 1; i <= 7; ++i) Vscanf("%ld", &mrx[i-1]); Vscanf("%*[^\n] "); Vscanf("%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))) { Vprintf("Allocation failure\n"); exit_status = -1; goto END; } for (i = 1; i <= nparx; ++i) Vscanf("%lf", &parx[i-1]); Vscanf("%*[^\n] "); /* Read model by which to filter series */ for (i = 1; i <= 3; ++i) Vscanf("%ld", &mr[i-1]); Vscanf("%*[^\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))) { Vprintf("Allocation failure\n"); exit_status = -1; goto END; } for (i = 1; i <= npar; ++i) Vscanf("%lf", &par[i-1]); Vscanf("%*[^\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)) ) { Vprintf("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)) ) { Vprintf("Allocation failure\n"); exit_status = -1; goto END; } g13byc(nser, &transfj, &fail); if (fail.code != NE_NOERROR) { Vprintf("Error from 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 g13bjc not to print parameters on entry */ options.list = FALSE; g13bjc(&arimaj, nser, &transfj, parx, npara, nx, iqxd, x, tdxxy, rms, mrxx, tdmrx, parxx, ldparx, tdparx, fva, fsd, &options, &fail); if (fail.code != NE_NOERROR) { Vprintf("Error from 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))) { Vprintf("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 g13bbc */ nmr = 10; nb = ny + MAX(mr[0] + mr[1], mr[2]); /* Allocate array b */ if ( !(b = NAG_ALLOC(nb, double)) ) { Vprintf("Allocation failure\n"); exit_status = -1; goto END; } /* Filter series by call to g13bbc */ g13bbc(y, ny, &transfv, &arimas, par, npar, cy, b, nb, &fail); if (fail.code != NE_NOERROR) { Vprintf("Error from g13bbc.\n%s\n", fail.message); exit_status = 1; goto END; } Vprintf(" Original Filtered\n"); Vprintf(" Backforecasts y-series series\n"); if (iqxd != 0) { ij = -iqxd; for (i = 1; i <= iqxd; ++i) { Vprintf("%8ld%17.1f%16.1f\n", ij, y[i-1], b[i-1]); ++ij; } Vprintf("\n"); Vprintf(" Filtered Filtered" " Filtered Filtered\n"); Vprintf(" series series" " series series\n"); for (i = iqxd + 1; i <= ny; i += 4) { for (ii = i; ii <= MIN(ny,i+3); ++ii) { Vprintf("%5ld", ii-iqxd); Vprintf("%10.1f", b[ii-1]); } Vprintf("\n"); } } } } END: /* Free the options structure used by g13bjc */ g13xzc(&options); if (b) NAG_FREE(b); if (fsd) NAG_FREE(fsd); if (fva) NAG_FREE(fva); if (par) NAG_FREE(par); if (parx) NAG_FREE(parx); if (x) NAG_FREE(x); if (y) NAG_FREE(y); if (rms) NAG_FREE(rms); if (parxx) NAG_FREE(parxx); if (mrxx) NAG_FREE(mrxx); return exit_status; }