/* nag_tsa_arma_filter (g13bac) Example Program. * * Copyright 2002 Numerical Algorithms Group. * * Mark 7, 2002. * Mark 7b revised, 2004. */ #include #include #include #include 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); 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; }