%============================================================
\subsection{\ML{} module installation}

Regular install.
\begin{mathematica}
Install["bj.ml"]
.
= ./bj.ml:1
\end{mathematica}

Install with IMSL's stdout redirected to stdin.
\begin{mathematica}
(
linkname = 9999;
While[
    Run["bj.ml -linkmode listen -linkname " <> ToString[linkname] <> " 1>&2 &"];
    Pause[1];
  Install[ToString[linkname], LinkMode->Connect] === $Failed, linkname -= 1];
  Links[]
)
.
= {LinkObject[9999@mecca, 1, 1]}
\end{mathematica}

Uninstall all.
\begin{mathematica}
Uninstall /@ Links[];
\end{mathematica}

%============================================================
\subsection{\texttt{GDATA}: Retrieve a commonly analyzed data set}
\begin{verbatim}
GDATA[idata_Integer, iprint:(0|1|2):0]
\end{verbatim}
\begin{mathematica}
x = #[[2]]& /@ (Take[#,{22,121+20}]& @ Transpose @ GDATA[2, 1]);
.
 The Wolfer Sunspot data.
 Reference:
    Anderson, T.W. (1971),
    The Statistical Analysis of Time Series,
    John Wiley and Sons, New York, page 660.
 This data set consists of 176 observations on 2 variables.


       First 10 rows of X
                1             2
  1        1749.0          80.9
  2        1750.0          83.4
  3        1751.0          47.7
  4        1752.0          47.8
  5        1753.0          30.7
  6        1754.0          12.2
  7        1755.0           9.6
  8        1756.0          10.2
  9        1757.0          32.4
 10        1758.0          47.6
\end{mathematica}

%============================================================
\subsection{\texttt{RNARM}:  Generate a
  time series from a specified ARMA model}
\begin{verbatim}
RNARM[nw_Integer, constant_?NumberQ, par_List, lagar_List, pma_List,
lagma_List, iadist:(0|1), avar_?NumberQ, a_List, wi_List].
Returns {w_List, a_List}.

Inputs:
NW - A Number of observations of the time series to generate.
     NW must be greater than or equal to one.
CONST - Overall constant.
PAR - Vector of length NPAR containing the autoregressive parameters.
LAGAR - Vector of length NPAR containing the order of the autoregressive
      parameters. The elements of LAGAR must be greater than or equal to one.
PMA - Vector of length NPMA containing the moving average parameters.
LAGMA - Vector of length NPMA containing the order of the moving average
      parameters. The elements of LAGMA must be greater than or equal to one.
IADIST - Option for normally distributed innovations.
  0: Innovations are generated from a normal distribution (white noise) with
     mean 0 and variance AVAR.
  1: Innovations are specifed by the user.
AVAR - Variance of the normal distribution, if used.  For IADIST = 0,
      AVAR is input; and for IADIST = 1, AVAR is unused.
A - Vector of length NW + max(LAGMA(.j.2)) containing the innovations.
      For IADIST = 1, A is input; and for IADIST = 0, A is output.
WI - Vector of length max(.LAGAR.(.i.>)) containing the initial values of
     the time series.

Outputs:
W - Vector of length NW containing the generated time series.
A - Vector of length NW + max(LAGMA(.j.2)) containing the innovations.
\end{verbatim}

\begin{mathematica}
(
{y, x2} = RNARM[ 1100, 0, {1.084, -0.2636}, {1, 2}, {700./1092} ,{1}, 0, 1.0, {}, {0., 0.}];
x = Drop[x2, 2];
)
\end{mathematica}

%============================================================
\subsection{\texttt{NSPE}: Compute preliminary
  estimates of the autoregressive and moving average parameters of an ARMA model}

\begin{mathematica}
(
imean = 0;
maxbc = 5;
lagar = Table[i, {i, 1, 2}];
lagma = Table[i, {i, 1, 2}];
lagma = {1};
npar = Length[lagar];
npma = Length[lagma];
)
\end{mathematica}

\begin{verbatim}
NSPE[w_List, imean:(0|1), wmean_?NumberQ, npar_Integer, npma_Integer,
  relerr_?NumberQ, maxit_Integer, iprint:(0|1):0].
Returns {wmean, constant, par_List, pma_List, avar}.
\end{verbatim}
\begin{mathematica}
{lwmean, lco, lpar, lpma, lavar} = NSPE[ x, imean, 0., npar, npma, 0.0, 0, 0]
.
= {0., 0., {3.94666, -0.000359466}, {0.253382}, 68922.8}
\end{mathematica}

%============================================================
\subsection{\texttt{NSLSE}: Compute least-squares estimates of parameters for a
nonseasonal ARMA model}
\begin{verbatim}
NSLSE[w_List, imean:(0|1), wmean_?NumberQ, par_List, lagar_List,
pma_List, lagma_List, maxbc_Integer, tolbc_?NumberQ, tolss_?NumberQ,
iprint:(0|1|2):0].
Returns {wmean, constant, par, pma, cov_?MatrixQ, a_List, avar}.
\end{verbatim}
\begin{mathematica}
{lwmean, lco, lpar, lpma, pcov, pa, lavar} =
  NSLSE[ x, imean, lwmean, lpar, lagar, lpma, lagma, maxbc, 0.0, .0, 2];
.
 ----------------------------------------------------------------------
 Iteration       1


        PAR
       1        2
  0.9420   0.0406

   PMA
  0.9912

 Residual SS (including backcasts) =  1085.305223333580
 Number of residuals               =    1099
 Number of backcasts               =       1
 ----------------------------------------------------------------------
 Final Results, Iteration       2

 Parameter             Estimate            Std. Error            t-ratio

    PAR
     1                 0.9420355            0.0323692           29.1028149
     2                 0.0405639            0.0289568            1.4008430

    PMA
     1                 0.9912233            0.0085595          115.8035382

 CONST =          0.0000000
 AVAR  =          0.9893393

 Residual SS (including backcasts) =       1085.3052233
 Number of residuals               =       1099

 Residual SS (excluding backcasts) =       1083.0391618
 Number of residuals               =       1098
\end{mathematica}

Calculate the t-ratios (the are not returned by NSLSE):
\begin{mathematica}
tratio = Prepend[ Join[ lpar, lpma], lwmean] /
  Prepend[ Sqrt[ MapIndexed[#1[[#2[[1]]]]&, pcov]], lavar]
.
= {0., 29.1028, 1.40084, 115.804}
\end{mathematica}

%============================================================
\subsection{\texttt{NSBJF}: Compute Box-Jenkins forecasts and their
associated probability limits for a nonseasonal ARMA model}
\begin{verbatim}
NSBJF[w_List, par_List, lagar_List, pma_List, lagma_List, iconst:(0|1),
constant_?NumberQ, avar_?NumberQ, alpha_?NumberQ, mxbkor_Integer,
mxlead_Integer, iprint:(0|1):0].
Returns fcst_?MatrixQ
\end{verbatim}

\begin{mathematica}
{out, problims, psi} = NSBJF[ x, lpar, lagar, lpma, lagma, 1, lco, lavar, 0.05, 0, Length[x], 0];
\end{mathematica}

%============================================================
\subsection{\texttt{LOFCF}: Perform lack-of-fit test for a univariate time series or transfer function given the appropriate correlation function}
\begin{verbatim}
LOFCF[nobs_Integer, lagmin_Integer, lagmax_Integer, cf_List, npfree_Integer].
Returns {Q, pvalue}.
\end{verbatim}
\begin{mathematica}
{npfree = npar + npma, lagmax = 25 + 0 Round[ Sqrt[N@Length[train]]], dof = lagmax - npfree}
.
= {3, 25, 22}
\end{mathematica}

\begin{mathematica}
residcor = ACF[pa, 1, imean, 0., lagmax, 0];
\end{mathematica}
\begin{mathematica}
ListPlot[residcor[[3]], PlotJoined->True];
\end{mathematica}

\begin{mathematica}
{Q, pvalue} = LOFCF[ Length[x], 1, lagmax, residcor[[3]], npfree]
.
= {30.8997, 0.901815}
\end{mathematica}

%============================================================
\subsection{\texttt{NSBJF}:  Compute Box-Jenkins forecasts and their
associated probability limits for a nonseasonal ARMA model}
\begin{verbatim}
NSBJF[w_List, par_List, lagar_List, pma_List, lagma_List, iconst:(0|1),
constant_?NumberQ, avar_?NumberQ, alpha_?NumberQ, mxbkor_Integer,
mxlead_Integer, iprint:(0|1):0].
Returns fcst_?MatrixQ
\end{verbatim}
\begin{mathematica}
{out, problims, psi} = NSBJF[ x, lpar, lagar, lpma, lagma, 1, lco, lavar, 0.05, 0, Length[x], 0];
\end{mathematica}


%============================================================
\subsection{\texttt{ARMME}: Compute
  method of moments estimates of the autoregressive parameters of an ARMA model}
\begin{verbatim}
ARMME[acv_List, npar_Integer, npma_Integer, iprint:(0|1):0].
Returns par_List
\end{verbatim}
\begin{mathematica}
amtest = ARMME[ artest[[2]], npar, npma, 1]
.
   Output PAR
      1       2
  3.947   0.000

= {3.94666, -0.000359466}
\end{mathematica}

%============================================================
\subsection{\texttt{MAMME}: Compute method of moments estimates of the moving
  average parameters of an ARMA model}
\begin{verbatim}
MAMME[acv_List, par_List, relerr_?NumberQ, maxit_Integer, npma_Integer,
iprint:(0|1):0].
Returns pma_List.
\end{verbatim}
\begin{mathematica}
MAMME[ artest[[2]], amtest, 0.0, 0, npma, 1]
.
 Output PMA from DMAMME/DM2MME
             0.2534

= {0.253382}
\end{mathematica}

%============================================================
\subsection{\texttt{RNSET}: Initialize a random seed for use in the IMSL random number generators}
\begin{verbatim}
RNSET[iseed_Integer]
\end{verbatim}
\begin{mathematica}
RNSET[123457]
\end{mathematica}


%============================================================
\subsection{\texttt{DIFF}: Difference a time series}
\begin{verbatim}
DIFF[z_List, iper_List, iord_List, iprint:(0|1):0].
\end{verbatim}
\begin{mathematica}
x = DIFF[ x, {1}, {3}];
\end{mathematica}

An equivalent function in \Mma.
\begin{mathematica}
diff[x_List, iper_List, iord_List] := Fold[ Function[{xi, p},
    Drop[ Nest[ RotateLeft[#, p[[1]] ] - #&, xi, p[[2]] ], -p[[1]] * p[[2]] ]],
  x, Transpose[{iper, iord}]]
\end{mathematica}


%============================================================
\subsection{\texttt{ACF}: Compute the sample autocorrelation function of a
stationary time series}
\begin{verbatim}
ACF[x_List, iseopt:(0|1|2), imean:(0|1), xmean_?NumberQ, maxlag_Integer,
iprint:(0|1|2|3):0].
Returns {xmean, acv_List, ac_List, seac_List}
\end{verbatim}

\begin{mathematica}
artest = ACF[ x, 1, imean, 0., 100, 1];
.
 Output from DACF/DA2F

 Mean     =    0.00000
 Variance =     4425.0
\end{mathematica}
\figi{sunspot}{Autocorrelation function of sunspot data.}

A rough sketch of an equivalent function in \Mma.
\begin{mathematica}
With[ {mu = Mean[x], n = N@Length[x], maxlag = 10, iseopt = 1},
  With[ {sigma = 1/n Table[ Sum[ (x[[t]] - mu) (x[[t+k]] - mu), {t, 1, n
      - k}], {k, 0, maxlag}]},
    With[ {rho = sigma / sigma[[1]]},
      hrho[i_Integer] := If[ Abs[i] <= maxlag, rho[[Abs[i]+1]], 0];
      With[ {std = Switch[iseopt,
               0, {},
               1, 1/n Table[ Sum[ hrho[i]^2 + hrho[i-k] hrho[i+k] -
                  4 hrho[i] hrho[k] hrho[i-k] + 2 hrho[i]^2 hrho[k]^2,
                  {i, -2 maxlag, 2 maxlag}], {k, 1, maxlag}],
               2, Table[ (n - k) / (n (n + 2)), {k, 1, maxlag}]]},
          {mu, sigma, rho, std}
]]]]
\end{mathematica}


%============================================================
\subsection{\texttt{PACF}: Compute the sample partial autocorrelation function of a
  stationary time series}
\begin{verbatim}
PACF[maxlag_Integer, ac_List].
Returns pac_List
\end{verbatim}

\begin{mathematica}
pacf = PACF[100, artest[[3]]];
\end{mathematica}


\vfill
\begin{rcslog}
$Log: bj.mma,v $
\Revision 1.5  1998/01/31 18:17:06  tjchol01
Before web.

\Revision 1.4  1997/10/10 21:00:02  tjchol01
Runs on mecca. Zurada exam.

\Revision 1.3  1996/09/30 16:01:27  tjchol01
Added |LOFCF|.

\Revision 1.2  1996/09/29 05:56:24  tjchol01
|bj.w| rewritten for ml.w.  Output redirection.

\Revision 1.1  1996/03/12 03:06:44  tjchol01
Initial revision
\end{rcslog}

% Local Variables:
% eval: (setenv "LM_LICENSE_FILE" "/usr/apps/vni/license/license.dat")
% End: