/* nag_rngs_varma_time_series (g05pcc) Example Program. * * Copyright 2001 Numerical Algorithms Group. * * Mark 7, 2001. */ #include #include #include #include int main(void) { /* Scalars */ Integer i, igen, ii, ip, iq, j, k, l, n, nr; Integer exit_status=0; NagError fail; Integer pdx, pdvar; Nag_OrderType order; /* Arrays */ double *phi=0, *r=0, *theta=0, *var=0, *x=0, *xmean=0; Integer iseed[4]; #ifdef NAG_COLUMN_MAJOR #define X(I,J) x[(J-1)*pdx + I - 1] #define VAR(I,J) var[(J-1)*pdvar + I - 1] order = Nag_ColMajor; #else #define X(I,J) x[(I-1)*pdx + J - 1] #define VAR(I,J) var[(I-1)*pdvar + J - 1] order = Nag_RowMajor; #endif INIT_FAIL(fail); Vprintf("nag_rngs_varma_time_series (g05pcc) Example Program Results\n\n"); /* Skip heading in data file */ Vscanf("%*[^\n] %ld%ld%ld%ld%*[^\n] ", &k, &ip, &iq, &n); nr = 600; /* Allocate memory */ if ( !(phi = NAG_ALLOC(k*k*ip, double)) || !(r = NAG_ALLOC(nr, double)) || !(theta = NAG_ALLOC(MAX(1,k*k*iq), double)) || !(var = NAG_ALLOC(k * k, double)) || !(x = NAG_ALLOC(k * n, double)) || !(xmean = NAG_ALLOC(k, double)) ) { Vprintf("Allocation failure\n"); exit_status = -1; goto END; } #ifdef NAG_COLUMN_MAJOR pdx = k; pdvar = k ; #else pdx = n; pdvar = k ; #endif if (n > 0 && n <= 100) { for (l = 0; l < ip; ++l) { for (i = 0; i < k; ++i) { ii = l * k * k + i; for (j = 0; j < k; ++j) { Vscanf("%lf", &phi[ii + k * j]); } Vscanf("%*[^\n] "); } } for (l = 0; l < iq; ++l) { for (i = 0; i < k; ++i) { ii = l * k * k + i; for (j = 0; j < k; ++j) Vscanf("%lf", &theta[ii + k * j]); Vscanf("%*[^\n] "); } } for (i = 0; i < k; ++i) { Vscanf("%lf", &xmean[i]); } Vscanf("%*[^\n] "); for (i = 1; i <= k; ++i) { for (j = 1; j <= i; ++j) Vscanf("%lf", &VAR(i,j)); Vscanf("%*[^\n] "); } /* Initialise the seed to a repeatable sequence */ iseed[0] = 1762543; iseed[1] = 9324783; iseed[2] = 4234401; iseed[3] = 742355; /* igen identifies the stream. */ igen = 1; /* nag_rngs_init_repeatable (g05kbc). * Initialize seeds of a given generator for random number * generating functions (that pass seeds explicitly) to give * a repeatable sequence */ nag_rngs_init_repeatable(&igen, iseed); /* nag_rngs_varma_time_series (g05pcc). * Generates a realisation of a multivariate time series * from a VARMA model */ nag_rngs_varma_time_series(order, 2, k, xmean, ip, phi, iq, theta, var, pdvar, n, x, pdx, igen, iseed, r, &fail); if (fail.code != NE_NOERROR) { Vprintf("Error from nag_rngs_varma_time_series (g05pcc).\n%s\n", fail.message); exit_status = 1; goto END; } Vprintf(" Realisation Number 1\n"); Vprintf("\n"); for (i = 1; i <= k; ++i) { Vprintf(" Series number %3ld\n", i); Vprintf(" -------------\n"); Vprintf("\n"); for (j = 1; j <= n; ++j) Vprintf("%8.3f%s", X(i,j), j%8 == 0 || j == n ?"\n":" "); Vprintf("\n"); } /* nag_rngs_varma_time_series (g05pcc), see above. */ nag_rngs_varma_time_series(order, 3, k, xmean, ip, phi, iq, theta, var, pdvar, n, x, pdx, igen, iseed, r, &fail); if (fail.code != NE_NOERROR) { Vprintf("Error from nag_rngs_varma_time_series (g05pcc).\n%s\n", fail.message); exit_status = 1; goto END; } Vprintf("\n\n"); Vprintf(" Realisation Number 2\n"); Vprintf("\n"); for (i = 1; i <= k; ++i) { Vprintf(" Series number %3ld\n", i); Vprintf(" -------------\n"); Vprintf("\n"); for (j = 1; j <= n; ++j) Vprintf("%8.3f%s", X(i,j), j%8 == 0 || j == n ?"\n":" "); Vprintf("\n"); } } END: if (phi) NAG_FREE(phi); if (r) NAG_FREE(r); if (theta) NAG_FREE(theta); if (var) NAG_FREE(var); if (x) NAG_FREE(x); if (xmean) NAG_FREE(xmean); return exit_status; }