Print[StyleForm[StringJoin[
        "Updated: October 30, 2003.",
        "Loaded UTA plus: ",
	"Statistics`DataManipulation`, ",
	"Statistics`DescriptiveStatistics`, ",
	"Statistics`NormalDistributions`. "
        ]
  ,FontSize->9,FontFamily->"Times", FontColor->
  RGBColor[0., 0.25,0.25]
   ]];

(* Author: A.I. McLeod *)
(* Copyright (c), A.I. McLeod, 2001 *)
(* UTA: Univariate Time Series Analysis Package *)
(* Summary:
A package for symbolic, numeric and graphic computations for univariate
time series analysis including arma models and spectral analysis.
This package requires auxilary packages: DatIO, QUtil.
*)


Needs["DatIO`", "d:/math/2003/DatIO.m"];
Needs["QUtil`", "d:/math/2003/QUtil.m"];

Needs["Statistics`DescriptiveStatistics`"];
Needs["Statistics`NormalDistribution`"];
Needs["Statistics`DataManipulation`"];

BeginPackage["UTA`"]

Needs["QUtil`", "d:/math/2003/QUtil.m"];
Needs["Statistics`DescriptiveStatistics`"];
Needs["Statistics`NormalDistribution`"];
Needs["Statistics`DataManipulation`"];
Needs["Statistics`MultinormalDistribution`"];
Needs["Statistics`ContinuousDistributions`"]; (* Panel..Fit *)
Needs["Statistics`Common`DistributionsCommon`"];
Needs["Graphics`MultipleListPlot`"];
Needs["Graphics`Animation`"];
Needs["LinearAlgebra`MatrixManipulation`"];

(*==================Initializations================================*)
Residuals;B;ResidualVariance;ARCoefficients;Pacf;
MeanZero;LoggedSdf;LoggedF;Normalized;
phi;MaximizedLogLikelihood;ChiSqNull;SignificanceLevel;
Determinant;
loglikelihoodar1;loglikelihoodrw;phimle;
SeedRandom[1];
Mu;Lambda;Omega;Delta;IVType;IVTime;IVErr;Loglikelihood;
WaldTest;GainTest;Wald;Gain
Step;Pulse;Ramp;Protect[Step];Protect[Pulse];Protect[Ramp];
\[Phi]hat;\[Mu]hat;avloglikelihoodar1;avloglikelihoodrw;\[Phi]opt;RLunitroot;

(*=========================Simulation==============================*)
SimulateTent::usage=
	"SimulateTent[n, z0] generates a sequence of length n from the tent map function 
with initial starting value z0 and z0 must lie in the open interval (0,1). 
The default value of z0 is 0.177221.";

SimulateGaussianIMA::usage=
"SimulateGaussianIMA[theta],n,Z1] simulates n consective observations
from an IMA(1) with unit innovation variance and initial value Z1.
The IMA(1) model equation is Z[t] = Z[t-1] + A[t] - theta*A[t-1],
where theta is the parameter and A[t] is white noise.";

SimulateAR1::usage=
"SimulateAR1[phi,n] simulates n successive observations from a stationary
Gaussian 
AR(1), z[t]=phi*z[t-1]+a[t], where innovation variance, var[a[t]]=1 and
Abs[phi]<1.";

ComputeMA::usage=
"ComputeMA[theta,a,n] is a compiled Mathematica function which 
computes a[[t]]-theta[[1]]*a[[t-1]]-...-theta[[q]]*a[[t-q]].";

SimulateGaussianMA::usage=
"SimulateGaussianMA[theta, n] simulates n consecutive observations from
a Gaussian MA(q) process
z[[t]] = a[[t]] - theta[[1]] a[[t-1]] - ... - theta[[q]] a[[t-q]],
where q=Length[theta] and a[[t]]'s are an NID(0,1) sequence.
In the MA(1) case, theta can be a number instead
of a vector. Not restricted to invertible processes.";

SimulateGaussianAR::usage=
"SimulateGaussianAR[phi, n] simulates n consecutive observations from
a stationary-causal Gaussian AR(p) process,
z[[t]] = a[[t]] + phi[[1]] z[[t-1]] - ... - phi[[p]] a[[t-p]],
where p=Length[phi] and a[[t]]'s are an NID(0,1) sequence.
In the AR(1) case, phi can be a number instead
of a vector.";

SimulateGaussianARMA::usage=
"SimulateGaussianARMA[phi, theta, n] simulates n consecutive observations from
a stationary-causal Gaussian ARMA(p,q) process, 
z[[t]] = a[[t]] + phi[[1]] z[[t-1]] - ... - phi[[p]] a[[t-p]]
 - theta[[1]] a[[t-1]] - ... - theta[[q]] a[[t-q]],
where p=Length[phi], q=Length[theta], a[[t]]'s are an NID(0,1) sequence and
phi and theta can be number vectors or scalars.
The MA component is not restricted to be invertible.";

SimulateGP::usage=
"SimulateGP[c] simulates a Gaussian process with autocovariance function c.
 The length of the realization is Length[c] which is assumed to be 
 greater than one.";

SimulateARIMA::usage = "SimulateARIMA[phi, theta, n] \
generates n consecutive random values from an ARMA process with parameters \
phi and theta and mean zero and unit innovation variance. 
SimulateSARIMA[phi, theta, sphi, theta, s] generates n consecutive 
random values from an SARIMA process with \
regular ARMA parameters phi and theta and seasonal ARMA parameters sphi, \
stheta and mean zero and unit innovation variance.
It takes options \
FractionalDifference and Differencing with both having default values 0. If \
FractionalDifference->d, Abs[d]<0.5 then fractionally differenced noise is \
generated. If Differencing->{d, zinit}, where d can take values 1, 2, ... and \
zinit is a vector of length d then an integrated process of order d is \
generated using the starting values in zinit.  Similarly if \
Differencing->{d,ds,s,zinit} where d,ds,s are positive integers and zinit is \
of length d+s*ds then an d-th order integrated and ds-th order \
seasonally-integrated process with seasonal period d is generated. If \
WarmUp>0, then n+WarmUp consecutive observations are generated and only the \
last n are returned.";

Options[SimulateARIMA] = {FractionalDifference -> 0, Differencing -> 0, 
     WarmUp -> 0};
 
Options[SimulateIA] = {IAType -> Step, Differencing -> 0, SigmaSquared -> 1};
 
SimulateIA /: SimulateIA::argserr = "All lists must be of length `1`.";
 
SimulateIA /: SimulateIA::usage = "SimulateIA[n, T, omega, delta, phi, theta] \
simulates an intervention which takes effect at time t=T in a series of \
length n.  The ARMA error noise parameters are phi={\!\(\[Phi]\_1\), ..., \
\!\(\[Phi]\_p\)} and theta={\!\(\[Theta]\_1\), ..., \!\(\[Theta]\_q\)} and \
the intervention parameters are omega = {\!\(\[Omega]\_0\), ..., \
\!\(\[Omega]\_g\)} and delta = {\!\(\[Delta]\_1\), ..., \!\(\[Delta]\_h\)}. \
The default is for a step intervention, IAType->Step but other possible types \
are IAType->Pulse and IAType->Ramp.  Here is a simple example to simulate a \
single pulse intervention with n=100, T=50 and \
\!\(\*Cell[TextData[Cell[BoxData[\n\(TraditionalForm\`\[Omega]\_0 = \
2\)]]]]\)and \!\(\*Cell[TextData[Cell[BoxData[\n\(TraditionalForm\`\[Phi]\_1 \
= 0.5\)]]]]\) SimulateIA[100,50,{2},{},{0.5},{},IAType->Pulse].  For \
convenience, in the case of a single intervention, if omega, delta, phi and \
theta are all vectors of length 1 or less then one can simply use the numeric \
values -- for example SimulateIA[100,50,1,0,0.5,0, IAType->Pulse] for the \
last example.  More generally multiple interventions may be simulated by \
making the parameter T a vector and the parameters omega, and delta a list of \
vectors.  By default all interventions are step.  To specify otherwise, in \
the multiple intervention case, you need specify a suitable vector of \
options.  Here is a simple example, involving simulate a series with two \
interventions the first a pulse and the second a step, \
SimulateIA[n,{50,60},{{2},{1}},{{0.7},{}},{0.5},{},IAType->{Pulse, Step}].  \
Integrated times series or models with differencing can be specified using \
the option Difference.  The default is no differencing, Differencing->0.  For \
regular differencing d times with d initial time series values in the vector \
zinit, set Differencing->{d, zinit}.  For regular differencing d times and \
seasonal differencing ds times with seasonal period s and d+ds*s initial time \
series values in the vector zinit, set Differencing->{d,ds,s,zinit}. The \
innovation variance is SigmaSquared which has default value 1.";

(*=========================SDF==============================*)
Options[CISdf]={Normalized->False};

OptimalPeriodogramSmoother::usage=
"OptimalPeriodogramSmoother[z] computes the optimal minimum risk Daniell smoother order.";

CISdf::usage=
"CISdf has two forms. CISdf[z] computes the 95% confidence intervals
for the spectral density function based on the optimal minimum risk
smoother.
CISdf[z,q] computes the 95% confidence interval for the spectral density 
function using an iterated Daniell smoother with half-window parameters given 
in q. The result is a matrix of 3 rows comprising: frequency, lower limit and 
upper limit respectively. 
CISdf takes one option with default: Normalized->False.
If Normalized->True, the normalized periodogram is used.";

SdfCheck::usage=
"Two forms of usage: SdfCheck[z, estsdf] or SdfCheck[z, estsdf, q].
SdfCheck[z, estsdf] computes the optimal minimum risk nonparametric
estimate of the sdf of data z to obtain 95% confidence intervals
which are then compared with the given sdf estimate provided by estsdf.
SdfCheck[z, estsdf, q] spectral analysis check for model of data z with given 
sdf, estsdf, using  
an iterated Daniell smoother with half-window widths in q. 
Normalized spectral density functions are used.
The spectral density estimate, estsdf,
can be obtained from the function 
SdfARMA[phi, theta, Normalized->True] 
after fitting a suitable model. 
SdfCheck takes options and defaults: LoggedF->False, 
LoggedSdf->False, PlotLabel->None.";

SdfPlot::usage=
"SdfPlot[z,q] plots estimated spectral density function. 
Options and defaults:
ARAlgorithm->\"Burg\",
Normalized->False,
ModelSelection->Automatic,
MeanZero->False,
Plot->True,
LoggedSdf -> False,
LoggedF->False, 
PlotRange->Automatic,
PlotLabel->None.
Two methods are 
available: nonparametric based on the iterated Daniell window and 
autoregressive spectral density estimation. The parameter q controls which 
method is used. If q is a vector of one or more integers, then the iterated 
Daniell smooth with half-width windows specified by q is used. For q={}, the 
raw periodogram is produced. If q is a single integer then the autoregressive 
spectral density estimator is used. In the autoregressive case, the default 
behaviour is that the best AR model of order between 0 and q is chosen using 
the AIC criterion and Burg parameter estimates. The other criterion available 
are the BIC or no selection. This is controlled by setting 
ModelSelection->BIC or ModelSelection->None. In the case of 
ModelSelection->None, an AR(q) is used for the spectral density estimate. The 
algorithm used is controlled by ARAlgorithm. ARAlgorithm->\"YuleWalker\" uses 
the Yule-Walker estimates instead of Burg estimates. 
If Plot->False, there is no plot and just the matrix with 2 columns comprising
the Fourier frequencies and the smoothed periodogram values is returned.
The axes scaling and
default values are: LoggedSdf->False, LoggedF->False. If LoggedSdf->True, a 
log base 2 scale is used for the vertical axis. If LoggedF->True, a log base 
2 scale is used for the horizontal scale. Related functions: Periodogram, 
SdfARMA, SdfARMAPlot, SdfARMAvsFARMAPlot.";

NPSdfPlot::usage=
"NPSdfPlot[z,q] nonparametric spectral density function estimate.
Options and defaults:
SdfSpan->Automatic,
Normalized->False,
MeanZero->False,
Plot->True,
LoggedSdf -> False,
LoggedF->False, 
PlotRange->Automatic,
PlotLabel->None.
If SdfSpan->Automatic the Daniell window if an automatic optimal
span is used. 
Otherwise if SdfSpan->q, where q is a vector,
of one or more integers, then the iterated 
Daniell smooth with half-width windows specified by q is used. 
For q={}, the raw periodogram is produced.
If Plot->False, there is no plot and just the matrix with 2 columns comprising
the Fourier frequencies and the smoothed periodogram values is returned.
The axes scaling and
default values are: LoggedSdf->False, LoggedF->False. If LoggedSdf->True, a 
log base 2 scale is used for the vertical axis. If LoggedF->True, a log base 
2 scale is used for the horizontal scale. Related functions: 
ARSdfPlot, Periodogram, 
SdfARMA, SdfARMAPlot, SdfARMAvsFARMAPlot, SdfCheck.";

ARSdfPlot::usage=
"ARSdfPlot[z,q] plots autoregressive spectral density function estimate.
Options and defaults:
EstimationAlgorithm->\"Burg\",
Normalized->False,
ModelSelection->Automatic,
MeanZero->False,
Plot->True,
LoggedSdf -> False,
LoggedF->False, 
PlotRange->Automatic,
PlotLabel->None.
ModelSelection->Automatic, uses the Bayes method for weighting all
AR(p) models for p=0,...,q.
If ModelSelection->\"AIC\" or ModelSelection->\"BIC\" then the best
AR(p), p=0,...,q is selected using the AIC or BIC criterion.
The particular estimation 
algorithm that is used is controlled by EstimationAlgorithm. 
The default is the Burg algorithm.
EstimationAlgorithm->\"YuleWalker\" uses 
the Yule-Walker estimates instead. 
If Plot->False, there is no plot and 
just the matrix with 2 columns comprising
the Fourier frequencies and the smoothed periodogram values is returned.
The axes scaling and
default values are: LoggedSdf->False, LoggedF->False. If LoggedSdf->True, a 
log base 2 scale is used for the vertical axis. If LoggedF->True, a log base 
2 scale is used for the horizontal scale. Related functions: 
NPSdfPlot, Periodogram, 
SdfARMA, SdfARMAPlot, SdfARMAvsFARMAPlot, SdfCheck.";

Options[SdfPlot]={ARAlgorithm->"Burg",
ModelSelection->Automatic,
MeanZero->False,
LoggedSdf -> False,
LoggedF->False, 
Normalized->False, 
Plot->True,
PlotRange->Automatic,
PlotLabel->None,
SamplingFrequency->1};

Options[ARSdfPlot]=
{EstimationAlgorithm->"Burg",
ModelSelection->Automatic,
MeanZero->False,LoggedSdf -> False,LoggedF->False, 
Normalized->False, Plot->True,
PlotRange->Automatic,
SamplingFrequency->1,
PlotLabel->None};

Options[NPSdfPlot]={
SdfSpan->Automatic,
MeanZero->False,LoggedSdf -> False,LoggedF->False, 
Normalized->False, Plot->True,
PlotRange->Automatic,
PlotLabel->None};

Options[SdfCheck]={LoggedSdf -> False, LoggedF->False, PlotLabel->None};

CumulativePeriodogramTest::usage=
"CumulativePeriodogramTest[z] plots the cumulative periodogram of time 
series z and 1% significance limits for testing for whiteness.";

NormalizedPeriodogram::usage=
"NormalizedPeriodogram[z] computes the normalized periodogram of 
time series z. Two rows are 
returned. The first row contains the frequencies and the second the 
normalized periodogram values. The periodogram is normalized by
dividing by the sample variance.";

Periodogram::usage=
"Periodogram[z] computes the periodogram of time series z. Two rows are 
returned. The first row contains the frequencies and the second the 
periodogram values.";

IteratedDaniellWeights::usage=
"IteratedDaniellWeights[q] for a given list of half-window width parameters 
for the modified Daniell smoother, the weights in the iterated smoother are 
returned.";

DaniellSmoother::usage=
"DaniellSmoother[p,q] smooths the periodogram in p with a modified Daniell 
smoother. 
The numeric vector p is the comprised of the periodogram values which
is the second element in the output from Periodogram, ie. p=Last[Periodogram[z]].
Note that when 
smoothing at frequencies close to 0 or 0.5 we use the symmetry inherent in 
the periodogram to wrap around so that frequency 1/n is used for -1/n and 
similarly 0.5+1/n is replaced by 0.5-1/n, etc., ie. periodic smoother.";


TacvfFDWN::usage=
"TacvfFDWN[d,MaxLag] computes the autocovariance function for fractionally 
differenced white noise with unit variance and parameter d at lags 0, 1, ..., 
MaxLag.";

SdfAR::usage = 
    "SdfAR[phi] evaluates the spectral density function of an 
autoregressive model, z[[t]] = a[[t]] + phi[[1]]*z[[t-1]] + ... + 
phi[[p]]*z[[t-p]], at frequencies 0,1/1024,...,511/1024 using the FFT. The 
parameter phi must be a number vector. The algorithm does not check if the 
parameters are in the admissible region. See InvertibleQ."; 

SdfMA::usage = 
    "SdfMA[theta] evaluates the spectral density function of a 
moving-average model, z[[t]] = a[[t]] - theta[[1]]*a[[t-1]] - ... - 
theta[[q]]*a[[t-q]], at frequencies 0,1/1024,...,511/1024  using the FFT. The 
parameter phi must be a number vector. The algorithm does not check if the 
parameters are in the admissible region. See InvertibleQ."; 

SdfARMA::usage = 
    "SdfARMA[phi,theta] evaluates the spectral density function of an 
autoregressive model, z[t] = z[[t]] = a[[t]] - theta[[1]]*a[[t-1]] - ... - 
theta[[q]]*a[[t-q]] + phi[[1]]*z[[t-1]] + ... + phi[[p]]*z[[t-p]], at 
frequencies m/1024,...,511/1024 using the FFT. Options and their defaults 
are: InnovationVariance->1,
BenchmarkTimeSeriesLength->1024, FractionalDifference->0,
Normalized->False. 
BenchmarkTimeSeriesLength is used to control the lowest frequency that
is shown. 
For a time series of length n, the lowest frequency in the periodogram
is 1/n.
Setting BenchmarkTimeSeriesLength->n means all frequencies shown in the
plot are greater than or equal to 1/n provided n<=1024.
The fast Fourier algorithm that is used to compute SdfARMA evaluates the
sdf at frequencies which are multiples of 1/1024.
By default, 
When the optional parameter 
BenchmarkTimeSeriesLength is set to some value less than 1024 then 
m=Ceiling[1024/BenchmarkTimeSeriesLength].  
When Normalized->True, the spectral density function is divided
by the theoretical variance so that the resulting function 
integrates to one.
The parameters phi and theta must 
be a number vectors. 
The algorithm does not check if the parameters are in 
the admissible region -- see InvertibleQ, ARMARedundancyQ and 
ARMAFactorsPlot. 
SdfARMA returns a matrix of the form 
{{\!\(f\_1\),\!\(sdf\_1\)},...}, where f is frequency and sdf is the spectral 
density at that frequency."; 

SdfFDWN::usage = 
    "SdfFDWN[d, n] evaluates the spectral density function for fractionally 
differenced white noise with parameter d at frequencies, 
Ceiling[1024/n]/1024, ..., 511/1024. The parameter n represents the minimum 
frequency that could be estimated by the periodogram in 
a series of length n.";

SdfARMAPlot::usage=
"SdfARMAPlot[phi,theta] plots the spectral density function of an 
autoregressive model, z[t] = z[[t]] = a[[t]] - theta[[1]]*a[[t-1]] - ... - 
theta[[q]]*a[[t-q]] ++ phi[[1]]*z[[t-1]] + ... + phi[[p]]*z[[t-p]]. Options 
and their defaults are: LoggedSdf->False, LoggedF->False, 
BenchmarkTimeSeriesLength->1024, FractionalDifference->0,
PlotRange->All, PlotLabel->None.
By default, BenchmarkTimeSeriesLength->1024 implies 
that m=1. If the optional parameter BenchmarkTimeSeriesLength is set to some 
value less than 1024 then m=Ceiling[1024/BenchmarkTimeSeriesLength]. If 
LoggedSdf->True, the log to base 2 of the spectral density is plotted. If 
LoggedF->True, the frequency axis is on a log base 2 scale. The parameters 
phi and theta must be a Number vectors. The algorithm does not check if the 
parameters are in the admissible region -- see InvertibleQ, ARMARedundancyQ 
and ARMAFactorsPlot.";

SdfARMAvsFARMAPlot::usage=
"SdfARMvsFARMAPlot[armaphi,armatheta,farmaphi,farmatheta,d] plots
the normalized spectral density functions from ARMA and 
fractional ARMA models together on the same graph. 
The normalized spectral density function is defined as the spectral density
function divided by the variance of the time series.
Options and their defaults 
are: LoggedSdf->False, LoggedF->False, BenchmarkTimeSeriesLength->1024,
PlotRange->All, PlotLabel->None, PLegend->{\"ARMA\",\"FARMA\"},
PStyle -> {RGBColor[1, 0, 1], RGBColor[0, 1, 0]}
By 
default, BenchmarkTimeSeriesLength->1024 implies that m=1. 
If the optional parameter 
BenchmarkTimeSeriesLength is set to some value less than 1024 then 
m=Ceiling[1024/BenchmarkTimeSeriesLength]. If LoggedSdf->True, the log to 
base 2 of the spectral density is plotted. If LoggedF->True, the frequency 
axis is on a log base 2 scale. The parameters phi and theta must be a Number 
vectors. 
The algorithm does not check if the parameters are in the admissible 
region -- see InvertibleQ, ARMARedundancyQ and ARMAFactorsPlot.";

Options[SdfARMAvsFARMAPlot] = 
{LoggedSdf -> False, 
PLegend->{ ARMA , FARMA },
LoggedF -> False, 
BenchmarkTimeSeriesLength -> 1024,
PlotRange->All, 
PlotLabel->None,
PStyle -> {RGBColor[1, 0, 1], RGBColor[0, 1, 0]}
};

Options[SdfARMA] = 
  {InnovationVariance->1, FractionalDifference -> 0, 
   BenchmarkTimeSeriesLength -> 1024, Normalized->False};
 
Options[SdfARMAPlot] = 
  {LoggedSdf -> False, LoggedF -> False, FractionalDifference -> 0, 
   BenchmarkTimeSeriesLength -> 1024, PlotRange->All, PlotLabel->None};

(*===========================Visualization==============================*)
BNPlot::usage=
"NBPlot[data] Produces a multipanel display of the normal 
probability plot and boxplot.";

BWPlot::usage=
	"BWPlot[data] box-and-whisker plot of data vector or data matrix.
In the case of a matrix the columns are taken as variables and the
plot is produced horizontally so the distributions of the variables
are most easily compared. The box goes from Q1 to Q3. The median
is shown as a point between Q1 and Q3. The whiskers extend out to
adjacent values which are the most extreme values in (Q1-r,Q3+r),
r=1.5(Q3-Q1). Data values outside this interval are plotted. Under
normality the probability of an outside value is about 0.7%.";

NPlot::usage=
	"NPlot[data] normal probability plot of data with quartile line.";

PhasePlot::usage=
	"PhasePlot[z, lag] produces a plot of z[t] vs z[t-lag] with successive 
values joined by a straight line. This plot is useful in visualizing the 
lagged dependence in the series.";

TimeSeriesPlot::usage=
	"TimeSeriesPlot[z,sublength,aspectratio] multipanel time series display. 
Each panel shows a subseries determined by the parameter sublength. 
The default for sublength is 200 and the
default for the aspect ratio is 0.25, somewhat less than 0.618 which is the 
standard Mathematica default. The aspect ratio should be chosen by the user 
to obtain the best visualization.";

AnimateTimeSeriesPlot::usage=
"AnimateTimeSeriesPlot[z,sublength,advance,aspectratio:0.25]
produces a sequence of time series plots which can be then
be run by an animation command (Ctrl-Y). The parameters are
sublength, the length of the time series shown on the plot;
advance, the number of observations skipped over between frames
which should be a minimum of 1 and a maximum of sublength-1 and
aspectratio, the ratio of the physical plot length divided by
the physical plot width.  The default value of aspectratio is 1.";

PlotARMAFactors::usage=
  "PlotARMAFactors[phi, theta] plots the reciprocal of the roots of the 
polynomials specified by 
1 - phi[[1]]*B - ... - phi[[Length[phi]]] B^Length[phi]
and
1 - theta[[1]]*B - ... - theta[[Length[theta]]] B^Length[theta]
in the complex plane along with the unit circle. The AR roots are
plotted as blue dots and the MA roots as brown dots. PlotARMAFactors[phi]
just plots the roots for the AR component.";

PlotSARMAFactors::usage=
  "PlotSARMAFactors[phi, theta, PHI, THETA, S] plots the reciprocal 
of the roots of the polynomials specified by 
1 - phi[[1]]*B - ... - phi[[Length[phi]]] B^Length[phi]
1 - theta[[1]]*B - ... - theta[[Length[theta]]] B^Length[theta] 
1 - PHI[[1]]*(B^S) - ... - PHI[[Length[PHI]]] (B^S)^Length[PHI]
and
1 - THETA[[1]]*(B^S) - ... - THETA[[Length[THETA]]] (B^S)^Length[THETA]
in the complex plane along with the unit circle. The AR roots and
MA roots are plotted as large blue dots and brown dots respectively. 
The seasonal AR roots and MA roots are plotted as small green dots 
and red dots.";


(*=========================TACF, etc.==============================*)
TacvfARMA::singular = 
    "Non-invertible AR or MA model or very close to it: the A Matrix is \
singular or ill-conditioned.";

Options[TacvfARMA] = {FractionalDifference -> 0, LagTruncation -> 512};

TacvfARMA::"unavail"="Symbolic fractional autocovariances are not available";

TacvfARMA::usage=
 "TacvfARMA[phi, theta, maxlag] evaluates the theoretical autocovariance
function at lags 0,1,...,maxlag for the ARMA model specified with
autoregressive component, phi, and moving-average component, theta.
The arguments phi and theta are vectors and maxlag is an
integer greater than or equal 0.
Two options are available to handle fractional arma models:
FractionalDifference -> 0 and LagTruncation -> 512.
When fractional difference is 0, we get the usual ARMA autocovariances
but for d, -0.5 < d < 0.5, fractional ARMA autocovariances are computed
using a high-order approximation which truncates the ARMA autocovariances
after the lag in indicated. Normally 512 is enough but it d is close to
0.5, an even larger value may be necessary to achieve reasonable accuracy.
";

CovarianceMatrixARMA::usage=
"CovarianceMatrixARMA[phi, theta, n] computes the n-by-n covariance matrix
of an ARMA model with parameters phi and theta. It takes the optional parameter
FractionalDifference.";

InverseCovarianceMatrixARMA::usage = "InverseCo\
varianceMatrixARMA[phi, theta, n] uses an exact efficient algorithm to \
determine the inverse of the covariance matrix of n successive observations \
from an ARMA process with parameters phi and theta. It takes the option \
Determinant with default: Determinant->False.  If Determinant->True, then a \
list is returned, the first element of which is the covariance determinant \
and the second is the inverse.";

SiddiquiMatrix::usage=
 "SiddiquiMatrix[phi, n] given phi a list of AR-parameters and n,
the n-by-n matrix which is the
exact inverse of the n-by-n covariance matrix of n successive
observations from the autoregression:
z[[t]] - phi[[1]]*z[[t-1] - ... - phi[[Length[phi]]] z[[t-Length[phi]]] = a[[t]],
where a[[t]] is white noise with unit variance is returned.
For AR(1) models, phi can be a real number instead of a list.";

Options[InverseCovarianceMatrixARMA] = {Determinant -> False};

Options[InformationMatrixARMA] := {FractionalDifference -> False};

InformationMatrixARMA::usage=
"InformationMatrixARMA[phi,theta] determines the Fisher information matrix
per observation for the parameters phi and theta in an ARMA model.
InformationMatrixARMA[phi,theta,phis,thetas,s] determines the Fisher information matrix
per observation for the parameters phi, theta, sphi and stheta in a seasonal ARMA model.
There is one option FractionalDifference. The default is FractionalDifference->False.
If the model contains a fractional differencing parameter then set FractionalDifference->True.
Note that the information matrix does not depend on the value of the fractional differencing
parameter.
";

InformationMatrixIFARMA::usage=
  "InformationMatrixIFARMA[phi,theta,d] determines the Fisher information matrix
for the parameters phi, theta and d in an IFARMA model.";

AcvfFD::usage =
"AcfFD[n, d] evaluates the acvf of a FARIMA(0,d,0) model
at lags 0,1,2,...,n using 13.2.8 & 13.2.9 of Brockwell and Davis.
Uses rule-based techique with auxiliary function gammaFD.";

(*=========================ACF==============================*)
AcfPlot::usage=
   "AcfPlot[r] plots the autocorrelation function, where
r is a vector of autocorrelations
at lags 1, ..., Length[r].
Any options are passed straight though to Graphics.";

VarianceRacfARMA::usage=
"VarianceRacfARMA[phi, theta, MaxLag, n] determines the covariance
matrix for the vector of residual autocorrelations at lags 
1, ... , MaxLag for the ARMA model specified by phi and theta which
is fit to a time series of length n. The variances may be obtained
by selecting the diagonal elements of this matrix using the
function Transpose[...,{1,1}].";

Options[Sacf] = {MaxLag -> 1, MeanZero -> False};

Sacf::usage=
"Sacf[z] returns the sample autocorrelations for time series z at
lags 1, ... , MaxLag.
It takes options with default values: MaxLag->1 and MeanZero->False. 
If MeanZero->False, then the sample mean
of z is subtracted and the sacf is computed for the corrected series.
If MeanZero->True, no sample mean correction is made.";

Sacvf::usage=
"Sacvf[z] returns the sample autocovariances for time series z at
lags 0, 1, ... , MaxLag.
It takes options MaxLag->1 and MeanZero->False. 
If MeanZero->False, then the sample mean
of z is subtracted and the sacvf is computed for the corrected series.";

Options[Sacvf] = {MaxLag -> 1, MeanZero -> False};

SacfPlot::usage=
"SacfPLot[z] plots the sample autocorrelation for time series z. 
It takes options and defaults: MaxLag->15, MeanZero->False, BartlettLargeLag->0, 
Bound->\"Individual\" and PlotLabel->None. The sacf is plotted 
for lags 1, ... , MaxLag. If MeanZero->False, then the sample mean
of z is subtracted and the sacf is computed for the corrected series.
If BartlettLargeLag->0, the approximate variance of 1/Length[z] is used
for the computation of the 95% benchmark limits. More generally when
BartlettLargeLag->Q, the series is assumed to be Q-dependent which 
means that data separated more than Q lags are uncorrelated and in this
case Bartlett's large lag formula is used to compute the 95% benchmark
limits. For Bound->\"Simultaneous\", 95% benchmark limits apply jointly to all
plotted autocorrelations.  If Bound->\"Individual\", the 95% benchmark limits are
for each individual autocorrelation. Otherwise, when Bound->None, no limits are shown.
If PlotLabel->\"Some text\" then the plot is labelled.";
 
Options[SacfPlot] = 
  {MeanZero -> False, MaxLag -> 15, BartlettLargeLag -> 0, 
   Bound -> "Individual",PlotLabel->None};

(*=========================PACF==============================*)
PacfDL::usage=
"PacfDL[c] computes the partial autocorrelations or partial regression 
coefficients for lags 1,2, ..., L given c a vector of the autocovariances at 
lags 0,1,...,L. Whether the partial autocorrelations or partial regression 
coefficients is produced is controlled by an option, LinearPredictor, which 
has a default setting of False. 
When LinearPredictor->False, the result is just the vector of partial autocorrelations.
When LinearPredictor->True, the result is list of rules for ARCoefficients, 
ResidualVariance, Pacf.
For improved efficiency, the algorithm uses compiled function in the case
of numeric autocovariances.";

Options[PacfDL]={LinearPredictor->False};

Options[PacfBayes]={LinearPredictor->False, MeanZero->False};

PacfBayes::usage=
"PacfBayes[z, p] computes the partial autocorrelations and/or 
fits an AR(p) model. A Bayesian weighted average of all AR
models up to order is used.
Whether the pacf or the fitted model is produced is controlled by
the option LinearPredictor which by default is False.";

PacfToAR::usage=
"PacfToAR[pi] converts the pacf in pi to coefficients in an
AR(p), p=Length[pi]";

ARToPacf::usage=
"ARToPacf[phi] converts AR(p) coefficients in phi to pacf.";

Options[PacfPlot] = 
  {Method->"Burg", MeanZero -> False, MaxLag -> 15, 
   SimultaneousBound -> True};

PacfPlot::usage=
"PacfPlot[z] plots the sample autocorrelation for time series z. 
It takes options and defaults: Method->\"Burg\", MaxLag->15, MeanZero->False, 
and
SimultaneousBound->True. The pacf is plotted 
for lags 1, ... , MaxLag. If MeanZero->False, then the sample mean
of z is subtracted and the pacf is computed for the mean corrected series. If 
SimultaneousBound->False, the 95% benchmark limits are
for each individual autocorrelation. Otherwise, when 
SimultaneousBound->True, the 95% bound applies to all plotted
autocorrelations. The Burg algorithm is the default. The Yule-Walker 
algorithm is specified with Method->\"YuleWalker\". 
Related Functions: PacfDL, PacfBurg, ICFromPacf.";

Options[PacfBurg] = 
  {MeanZero -> False, LinearPredictor -> False, ModelSelection -> False};
 
Options[ICFromPacf] = {Criterion -> "AIC", Table->False,Plot -> False};

PacfBurg::usage=
"PacfBurg[z,MaxLag] computes the partial autocorrelations at lags 1, ..., 
MaxLag using the Burg algorithm. Options and defaults: MeanZero->False, 
LinearPredictor->False. If LinearPredictor->True, an AR(MaxLag) is fit and 
the result is list of rules for ARCoefficients, ResidualVariance, Pacf.";

ICFromPacf::usage=
"ICFromPacf[n,pacf,opts] given a pacf sequence and n, the length of the time
series, the optimal minimum AIC AR model order is returned. The AIC is the 
default, Criterion->\"AIC\". Other values allowed are \"BIC\". The other 
options and their defaults: Plot->False, Table->False. 
If Plot->True, the a plot of the IC 
criterion as a function of the model order k is produced.
If Table->True a vector is returned contain the IC values for
k=0,...,Length[pacf] and the option Plot is ignored.
";

(*=========================MA(1) Estimation==============================*)
Options[MA1Fit] = {MeanZero -> False};
 
MA1Fit::usage = 
"MA1Fit[z] returns the exact MLE for \[Theta] in 
the MA(1) model, z[t] = a[t] - \[Theta] a[t-1], where a[t] is Gaussian white 
noise.  It takes option MeanZero with default MeanZero->False. If 
MeanZero->False, the sample mean is subtracted from z and then the exact 
likelihood equation is solved for a mean-zero process. Otherwise if 
MeanZero->True, no sample mean correction is made.";
 
(*=========================AR(1) Estimation==============================*)

Options[MarkovModel] = {MeanZero -> False, ErrorTolerance -> 1/100000, 
     MeanAverage -> False};

MarkovModel::usage=
"MarkovModel[z] fits exact mle for mean and ar(1) coefficient.
Returns a list of rules for estimates and the average loglikelihoods
of the optimized ar(1) model and the random walk model; 
the optimal mle based on the average loglikelihood and
RL, the relative likelihood based of unit root. Options and defaults
are MeanZero->True, ErrorTolerance->1/100000, MeanAverage->False. 
This algorithm uses compiled Mathematica functions for maximum
computational efficiency.";

AR1Fit::usage = 
"AR1Fit[z] exact Gaussian maximum likelihood estimate for the parameter 
in a stationary AR(1) for time series z obtaining by solving the cubic
likelihood equation. It takes option MeanZero with default MeanZero->False.
If MeanZero->True, the sample mean is subtracted from z and then the 
exact likelihood equation is solved for a mean-zero process.";

Options[AR1Fit] = {MeanZero -> False};

(*=========================AR Estimation==============================*)
OrdinaryLeastSquaresARFit::usage=
"OrdinaryLeastSquareARFit[z, p] fits an AR(p) to the time series in vector z. Ordinary
least squares is used. The parameter vector is returned.";

ResidualsAR::usage=
"ResidualsAR[z, phihat] given the time series z[[1]], ... , z[[n]] and
the parameters phihat[[1]], ..., phihat[[p]], where n=Length[z] and
p=Length[phihat], the n-p residuals, 
a[[t]]=z[[t]]-phihat[[1]]*z[[t-1]]-...-phihat[[p]]*z[[t-p]],
t=p+1, ..., n, are calculated.";

DetAR::usage=
"DetAR[phi] determinant of an AR process with parameter vector phi.";
 
Options[LoglikelihoodAR] = {MeanZero -> False};

LoglikelihoodAR::usage = "LoglikelihoodAR[phi,z] computes 
the exact loglikelihood for the AR(p) with parameters phi and observed data z.
It takes option MeanZero with default value MeanZero->False. 
If MeanZero->False, then the sample mean of z is subtracted.
Otherwise if MeanZero->True, the mean of z is assumed to be zero.
The loglikelihood returned is -1/2*Log[g]-n/2 Log[S/n], where n=Length[z] 
g is the covariance determinant and S
is the sum of the standardized prediction errors.";

Options[ChampernowneD] = {MeanZero -> False};
 
ChampernowneD::usage = "ChampernowneD[z,p] returns the
(p+1)-by-(p+1) matrix with (i,j) entry
z\[LeftDoubleBracket]i\[RightDoubleBracket]*z\[LeftDoubleBracket]j\[RightDoub\
leBracket] + ... + \
z\[LeftDoubleBracket]n-j\[RightDoubleBracket]*z\[LeftDoubleBracket]n-i\[Right\
DoubleBracket].  It takes option MeanZero->False. If MeanZero->False, then \
the sample mean of z is subtracted.";

Options[YuleWalkerEstimates]={MeanZero->False};

YuleWalkerEstimates::usage=
	"YuleWalkerEstimates[z, p] computes the Yule-Walker estimates 
of the autoregressive 
parameters in an AR(p) model. The optional parameter and its default is 
MeanZero->False.";

Options[ChampernowneLS] = {MeanZero -> False, Lags -> Automatic, 
     Constrained -> False};
 
ChampernowneLS/: ChampernowneLS::inadmissible = 
     "Inadmissible AR estimates. Shrunken admissible result produced.";
 
ChampernowneLS/: ChampernowneLS::usage = "ChampernowneLS[z,p] computes the \
Champernowne LS estimates for the AR(p). It takes option MeanZero->False, \
Lags->Automatic, Constrained->False. If MeanZero->False, then the sample mean \
of z is subtracted. If Lags->{\!\(p\_1\),...,\!\(p\_k\)} then a subset \
AR(\!\(p\_k\)) is fitted using only lags \!\(p\_1\) ..., \!\(p\_k\). If \
Constrained->True then if inadmissible AR estimates are obtained the result \
is shrunk towards the origin until an admissible result is obtained.";

(*=========================ARMA Estimatoin==============================*)
DetARMA::usage=
"DetARMA[phi, theta] returns the asymptotic covariance determinant
of the standardized ARMA time series Z[1]/s, Z[2]/s, ... where
s is the standard deviation.";

Options[LoglikelihoodGP] = {MeanZero -> False};

LoglikelihoodGP::usage=
"LoglikelihoodGP[c, z] 
computes the exact loglikelihood for a time series with 
data z and autocovariance function c. 
If c has length less than z, it is padded with zeros. 
Both c and x are assumed to be nonempty numeric vectors.
It takes option MeanZero with default value MeanZero->False. 
If MeanZero->False, then the sample mean of z is subtracted.
Otherwise if MeanZero->True, the mean of z is assumed to be zero.
The loglikelihood returned is -1/2*Log[g]-n/2 Log[S/n], where n=Length[z] 
g is the covariance determinant and S
is the sum of the standardized prediction errors.";

ResidualsGP::usage=
"ResidualsGP[c, z] computes the standardized residuals for a GLP. It
takes the same options as LoglikelihoodGP.";

Options[LoglikelihoodARMA] = {MeanZero -> False, FractionalDifference->0};
 
LoglikelihoodARMA::usage = "LoglikelihoodARMA[phi,theta,z\
] evaluates the concentrated loglikelihood of an ARMA model with parameters \
phi and theta and data z.  It takes the option, MeanZero.  For the default \
MeanZero->False, the data is corrected by the sample mean. When \
MeanZero->True, no correction is made.";
 
(*=========================Panel Time Series==============================*)
PanelTimeSeriesQ::usage=
"PanelTimeSeriesQ[z] tests if z is a list each component of which is a vector.";
SARToAR::usage = "SARToAR[phi,sphi, s] returns gphi, where gphi is \
the ar parameter in the regular AR.";

FitPanelTimeSeriesAR1::usage=
"FitPanelTimeSeriesAR1[z] computes MLE and LR-test for panel time series.";

FitPanelTimeSeriesAR1::constant=
"Constant series in component deleted: `1` components remaining.";

LoglikelihoodPanelTimeSeries::usage=
"LoglikelihoodPanelTimeSeries[c, z] 
computes the exact loglikelihood for a panel time series.
The theoretical autocovariance function for the longest component
panel is specified by a numeric vector c.
The panel time series data is specified by z which is a list such that each
component of the list corresponds to one of the time series panels.";


(*=========================Intervention Analysis==========================*)

Options[IADesign] = {IAType -> Step, Differencing -> 0, Delta1 -> 0};
 
IADesign /: IADesign::invalidia = "Invalid IAType-> `1`, reset to Step.";
 
IADesign /: IADesign::notAvail = "Option for delta operator is not available \
for multiple interventions. Ignored.";
 
IADesign /: IADesign::usage = "IADesign[n,T,g] returns the design matrix, of \
dimensions n-by-g corresponding to an intervention in a series of length n \
which takes effect at time T and has an omega-operator term of degree \
g\[GreaterEqual]0.  The default, IAType->Step is for step interventions but \
pulse and ramp interventions can be used with IAType->Pulse and IAType->Ramp. \
 Noise terms with differencing are also supported. The default is with no \
differencing parameter, Differencing->0.  For regular differencing use, \
Differencing->d, where d=1,2 etc. For the seasonal case, use \
Differencing->{d, ds, s} where d specifies the regular differencing and ds \
specifies the degree of seasonal differencing and s is the seasonal period. \
In addition, IADesign, can be used with a known delta-operator, \
IADesign[n,T,g,Delta1->\[Delta]] returns the design matrix corresponding to a \
known delta operator of the form (1-\[Delta] B), where Abs[\[Delta]]<1.  The \
default value for Delta1 is, Delta1->0.";

Options[IAFit1] = {Verbose -> True, Noise -> "AR1", IAType -> Step, 
     Differencing -> 0, Delta1 -> 0, MaxitIAFit1 -> 100, 
     ErrorTolerance -> 1/10000};
 
IAFit1 /: IAFit1::maxiter = "Warn:Maximum number of iterations reached: `1` \
iterations. Abs rel err = `2`";
 
IAFit1 /: IAFit1::usage = "IAFit1[z,T,g] fits a single intervention analysis \
model with starting time, t=T, and omega operator of degree g to the series \
z.  The default intervention is IAType->Step, a step intervention but other \
options include IAType->Pulse and IAType->Ramp.  The default noise model is \
Noise->\"AR1\", with options Noise->\"NID\" and Noise->\"MA1\".  Additional \
default options are MeanZero->True and Verbose->True.  You can set \
Verbose->False, to stop printout of the iterations or Verbose->m to printout \
every m-th iteration.  Printouts are in blue font. Models with known \
\!\(\[Delta]\_1\) can be fit by specifying Delta1->\!\(\[Delta]\_1\).  The \
error tolerance has default, ErrorTolerance->\!\(10\^\(-4\)\), which is the \
maximum relative error in the estimation of the parameters.  MaxitIAFit1 \
iterations are done, default is MaxitIAFit1->100.";

Options[IAFit] = {VerboseBrent -> False, VerboseIAFit1 -> False, 
     Noise -> "AR1", IAType -> Step, Differencing -> 0, 
     EstimateDelta -> False};
 
IAFit /: IAFit::usage = "IAFit[z,T,g] fits a single intervention analysis \
model with starting time, t=T, and omega operator of degree g to the series \
z.  The default intervention is IAType->Step, a step intervention but other \
options include IAType->Pulse and IAType->Ramp.  The default noise model is \
Noise->\"AR1\", with options Noise->\"IID\" and Noise->\"MA1\".  Additional \
options and their default values are  VerboseBrent->False, \
VerboseIAFit1->Falsle, Differencing->0 and Delta->False.  Set Differencing->d \
for regular differencing and Differencing->{d,ds,s} for more general regular \
and/or seasonal differencing.  Set EstimateDelta->True to estimate the delta \
parameter.  VerboseBrent->True, to see the progress in the Brent optimization \
(magenta font).  VerboseIAFit1->False is default but to see the progress in \
IAFit1 iterations, set VerboseIAFit1->True or VerboseIAFit1->m to see every \
iteration or every m-th iteration (blue font).";

Options[TestIA] = {Test -> Both, IAType -> Step, TestType -> TwoSided, 
     SigmaSquared -> 1, PierceApproximation -> False};
 
TestIA /: TestIA::singular = "Information matrix is singular.";
 
TestIA /: TestIA::usage = "TestIA[n,T,omega,delta,phi,theta] compute the \
attained significance level, n=length of series, T=time of intervention, \
omega = numerator in TF, delta = denominator in TF, phi = ar, theta = ma.  \
Options and defaults are Test->Both (options: WaldTest, GainTest), \
Level->0.05, IAType->Step, TestType->TwoSided, SigmaSquared->1.";

PowerIA /: PowerIA::singular = "Information matrix is singular.";

Options[PowerIA] = {Test -> Both, Level -> 0.05, IAType -> Step, 
     TestType -> TwoSided, SigmaSquared -> 1, PierceApproximation -> False};
 
PowerIA /: PowerIA::usage = "PowerIA[n,T,omega,delta,phi,theta] compute the \
power in IA with parameters, n=length of series, T=time of intervention, \
omega = numerator in TF, delta = denominator in TF, phi = ar, theta = ma.  \
Options and defaults are Test->Both (options: WaldTest, GainTest), \
Level->0.05, IAType->Step, TestType->TwoSided, SigmaSquared->1.";
 
Options[ICalc] = {IAType -> Step};
 
ICalc /: ICalc::usage = "ICalc[n,T,omega,delta] computes the output of an \
intervention which takes effect at time t=T in a series of length n.  The \
intervention parameters are omega = {\!\(\[Omega]\_0\), ..., \
\!\(\[Omega]\_g\)} and delta = {\!\(\[Delta]\_1\), ..., \!\(\[Delta]\_h\)}. \
The default is for a step intervention, IAType->Step but other possible \
types are IAType->Pulse and IAType->Ramp.";

TFCalc /: TFCalc::usage = "TFCalc[omega,delta,x] computes the output of a \
transfer function defined by z[t] = (omega(B)/delta(B))x[t], where omega(B) \
and delta(B) are operators defined by \
omega(B)=omega\[LeftDoubleBracket]0\[RightDoubleBracket] + ... + \
omega\[LeftDoubleBracket]g\[RightDoubleBracket]B^g and delta(B) = 1 - \
delta\[LeftDoubleBracket]1\[RightDoubleBracket]B - ... - \
delta\[LeftDoubleBracket]h\[RightDoubleBracket]B^h, B is the backshift \
operator on t.  The input parameters are vectors \
omega={omega\[LeftDoubleBracket]0\[RightDoubleBracket], ..., \
omega\[LeftDoubleBracket]gp1]}, \
delta={delta\[LeftDoubleBracket]1\[RightDoubleBracket], ..., \
delta\[LeftDoubleBracket]h\[RightDoubleBracket]} and \
x={x\[LeftDoubleBracket]1\[RightDoubleBracket],...,x\[LeftDoubleBracket]n\[Ri\
ghtDoubleBracket]}.  The initial values are handled by setting x[t] = x[1] \
for all t<1.";

Options[InfMatIAP] = {SigmaSquared -> 1, X -> {}, 
        PierceApproximation -> False, IAType -> Step};

Options[InformationMatrixIA] = {SigmaSquared -> 1, X -> {}, 
     PierceApproximation -> False, IAType -> Step};
 
InformationMatrixIA /: InformationMatrixIA::usage = "InformationMatrixIA[n, \
T, omega, delta, phi, theta] computes the Fisher information matrix for the \
intervention analysis model with n successive observations, T is the \
intervention time point, omega is a vector of length g+1 containing the \
coefficients (omega[0], omega[1], ..., omega[g]) and delta is of length h \
containing (delta[1], ..., delta[h]), phi=(phi[1],...,phi[p]) and  \
theta=(theta[1],...,theta[q]). The intervention model is defined by z[t] = mu \
+ omega(B)/delta(B) I(t,T) + theta(B)/phi(B) a[t], where a[t] is white noise \
with unit variance, B is the backshift operator and the operators omega(B), \
delta(B), phi(B) and theta(B) are defined by omega(B) = omega[0] + omega[1]B \
+...+ omega[g]B^g, delta(B) = delta[0] + delta[1]B +...+ delta[g]B^h, phi(B) \
= phi[0] + phi[1]B +...+ phi[g]B^p, theta(B) = theta[0] + theta[1]B +...+ \
theta[g]B^q. The order of parameters in the information matrix is (mu,  \
omega[0], omega[1], ..., omega[g], delta[1], ..., delta[g]).  For the case of \
multiple interventions, the parameters omega and delta are lists of vectors.  \
This function, InformationMatrixIA has five optional arguments.  See \
Options[InformationMatrixIA] for settings.  One optional argument generalizes \
to the case where a[t] has variance \
\!\(\*Cell[TextData[Cell[BoxData[\n\(TraditionalForm\`\[Sigma]\^2\)]]]]\) \
another optional argument allows for linear covariates, X, a third allows for \
ARIMA and seasonal ARIMA models with differencing, a fourth option is for the \
mean-zero case and the fifth option is to use the Pierce approximation.  The \
default values for these optional arguments are SigmaSquared\[Rule]1,  \
X\[Rule]{} and Differencing->0.  More generally you can set \
SigmaSquared\[Rule]\!\(\[Sigma]\^2\) and X\[Rule]to a suitable n by M matrix \
of covariates.    If regular differencing to order d is included then set \
Differencing\[Rule]d and if seasonal differencing is present then use \
Differencing\[Rule]{d,ds,s}, where ds indicates the order of seasonal \
differencing and s is the seasonal period.  The default is MeanZero->False \
which can be reset to MeanZero->True if appropriate.  Finally, \
PierceApproximation->False can be set to PierceApproximation->True.";

SEGain /: SEGain::usage = "SEGain[omega,delta,imatrix] computes the standard \
error of the gain coefficient in an IA model with transfer-function \
parameters omega and delta and Fisher information matrix imatrix.  Note \
imatrix is obtained from InformationMatrixIA.";

VarianceGain /: VarianceGain::usage = "VarianceGain[omega,delta,imatrix] \
computes the variance of the gain coefficient in an IA model with \
transfer-function parameters omega and delta and Fisher information matrix \
imatrix.  Note imatrix is obtained from InformationMatrixIA.";

(*=========================INPUT/OUTPUT==============================*)

ReadARMAParameters::usage=
   "ReadARMAParameters[file] reads in the vector
of parameters {{phi_1,...,phi_p},{theta_1,...,theta_q}} from
a file whose first two values are p and q and whose remaining
values are phi_1,...,phi_p,theta_1,...,theta_q.";

ReadSARMAParameters::usage=
   "ReadSARMAParameters[file] reads in the vector
of parameters {{phi_1,...,phi_p},{theta_1,...,theta_q},{PHI_1,...,PHI_P},
{THETA_1,...,THETA_Q}} from a file whose first five values are p,q,P,Q and 
S and whose remaining values are phi_1,...,phi_p,theta_1,...,theta_q,
PHI_1,...,PHI_P,THETA_1,...,THETA_Q.";

(*=========================Data Transformation==============================*)
PowerTransformation::usage=
	"PowerTransformation[z, lambda] data z are transformed according to lambda. 
If lambda is 1, no transformation. If lambda is zero, logs are taken. 
Otherwise, z^lambda.";

TSIntegrate::usage = "TSIntegrate[w,d,ds,s,zinit] integrates 
the time series w using initial values zinit and with parameters d, ds, s. 
TSIntegrate[w,d,zinit] and TSIntegrate[z,ds,s,zinit] may be used for purely
regular or seasonal integration respectively.
TSIntegrate reverses TSDifference.";

TSDifference::usage=
"TSDifference[z,d,s] differences or seasonally differences a time series d 
times. The default value of d is 1. 
The default value of s, the seasonal period, is 1.
TSDifference[z,d,ds,s] difference d times followed by seasonal differencing
ds times with seasonal period s. Default values d:1, ds:0, s:12.";

(*=========================ARMA Utilities==============================*)
InvertibleQ::usage=
  "InvertibleQ[phi] returns True or False according as the polynomial specified
by 1 - phi[[1]]*B - ... - phi[[Length[phi]]] B^Length[phi] has all roots
outside the unit circle.";

SARMAToARMA::usage = "SARMAToARMA[phi,theta,sphi,stheta,s] returns gphi 
and gtheta, where gphi and gtheta are the 
parameters in the regular ARMA.";

ARRoots::usage=
  "ARRoots[phi] finds the roots of the polynomial specified
by 1 - phi[[1]]*B - ... - phi[[Length[phi]]] B^Length[phi].";

GenerateARMAPolynomial::usage=
  "GenerateARMAPolynomial[x] returns coefficients of a polynomial whose roots
 are the reciprocals of the elements in the list x.";

TestARMARedundancy::usage=
  "TestARMARedundancy[phi,theta] determines if the ARMA model specified by
the parameters phi and theta is redundant. If it is redundant, the number
of common roots is noted. Euclid's algorithm is used.";

AuxiliaryMatrixARMA::usage=
  "AuxiliaryMatrixARMA[phi,theta] determines the auxiliary matrix of the
ARMA model specified by phi and theta.";

ARMARedundancy::usage=
 "ARMARedundancy[phi, theta] given vectors of the AR and MA coefficients,
phi and theta respectively, the equation which gives a necessary
and sufficient condition for model redundancy is returned.";

Toeplitz::usage=
  "Toeplitz[x] forms the Toeplitz matrix from vector x.";

ARPolynomial::usage=
  "ARPolynomial[phi, x] = 1 - phi[1] x - ... - phi[p] x^p, where
p = Length[phi].";

Options[ImpulseCoefficientsARMA]={Difference->0, SeasonalPeriod->1,
SeasonalDifference->0};

Options[InvertedCoefficientsARMA]={Difference->0, SeasonalPeriod->1,
SeasonalDifference->0};

ToPhiStar::usage=
"ToPhiStar[phi, d, ds, s] coefficients of the generalized autoregressive 
component obtained by multiplying the regular component times the 
differencing and seasonal differencing components.";

ExpandToSeasonal::usage=
"ExpandToSeasonal[phis, s] coefficients in expanded seasonal operator.";

ARInvert::usage=
"ARInvert[phi, MaxLag] given a list phi of the AR parameters, the 
first m coefficients in the MA approximation are calculated.  Specifically	
the vector {psi[1], ..., psi[MaxLag]}, where
1/(1-phi[[1]]B-phi[[2]]B^2-...-phi[[p]]B^Length[phi]) = 
1+psi[1]B+...+psi[[m]]B^m + ... .";

ImpulseCoefficientsARMA::usage=
"ImpulseCoefficientsARMA[phi, theta, MaxLag] returns the vector 
{1,psi[[1]],psi[[2]],...,psi[[MaxLag]]} where the psi's are the coefficients 
are the a's in the expansion, 
z[[t]]=a[[t]]+psi[[1]]*a[[t-1]]+psi[[2]]*a[[t-2]]+... . Takes options and 
defaults, Difference->0, SeasonalPeriod->0, SeasonalDifference->0.";

InvertedCoefficientsARMA::usage=
"InvertedCoefficientsARMA[phi, theta, MaxLag] returns the vector 
{1,pi[[1]],pi[[2]],...,pi[[MaxLag]]} where the pi's are the coefficients are 
the z's in the expansion, z[[t]]=a[[t]]+pi[[1]]*z[[t-1]]+pi[[2]]*z[[t-2]]+... 
. Takes options and defaults, Difference->0, SeasonalPeriod->0, 
SeasonalDifference->0.";

FDtoGLP::usage = 
"FDtoGLP[d, L] returns the coefficients \!\(\[Psi]\_1\) , ... ,\!\(
\[Psi]\_L\) for an approximation \!\(z\_t\) = 
\!\(a\_t\)+\!\(\[Psi]\_1\)\!\(a
\_\(t - 1\)\)+...\!\(\[Psi]\_L\)\!\(a\_\(t - L\)\) to fractionally 
differenced white noise with parameter d."; 

FDtoAR::usage = 
"FDtoAR[d, L] returns the coefficients \!\(\[Phi]\_1\) , ... ,\!\(
\[Phi]\_L\) in an AR(L) approximation to fractionally differenced white noise 
with parameter d.";


(*=========================Utilities==============================*)
InformationMatrixToSE /: InformationMatrixToSE::usage = "InformationMatrixToS\
E[imatrix] determines the standard errors corresponding to the given Fisher \
information matrix.";

InformationMatrixToCorrelationMatrix /: 
    InformationMatrixToCorrelationMatrix::usage = "InformationMatrixToCorrela\
tionMatrix[imatrix] determines the correlation matrix corresponding to the \
given Fisher information matrix.";
 
IndexOfMinimum::usage=
"IndexOfMinimum[x] finds the index position of the minimum value in a vector 
of numbers x. In the case of ties, only the first minimum index is returned.";

Average::usage=
"Average[x] average of all elements in x.";
 
Options[Brent] = {ErrorTolerance -> 1/100000, GridSearch -> 10, 
     MaxIterBrent -> 40, Verbose -> False};
 
Brent /: Brent::Errfunc = 
     "Error: Objective function must have one argument, `1`";
 
Brent /: Brent::imaginary = 
     "Error: Objective function contains symbolic or imaginary component.";
 
Brent /: Brent::maxiter = "Warn: Maximum number of iterations reached: `1` \
iterations. Abs rel err = `2`";
 
Brent::usage = "Brent[func, xmin, xmax] evaluate argmax func[x] on \
the interval [xmin, xmax] using a robust combination of \ 
grid search, inverse quadratic interpolation and Fibonnaci search. \
Brent takes three \
optimal arguments and their default values are: \
ErrorTolerance->\!\(10\^\(-5\)\), GridSearch->10, MaxIterBrent->25 and \
Verbose->False. The  value of ErrorTolerance causes the algorithm to \
terminate when the relative width of the bracketing interval is less than \
ErrorTolerance.  If we set GridSearch->n, the initial bracketing interval \
[xmin, xmax] is initially divided into a grid on n equi-distant points and \
the best available bracketing interval is used to start the algorithm. \
MaxIterBrent sets the maximum number of iterations of the inverse quadratic \
interpolation algorithm.  Verbose->True causes a detailed printout of the \
progress -- magneta is used to indicate more clearly the output produced by \
Brent since the func itself may produce other output.";

(********************************************************************) 
(********************************************************************) 

Begin["`Private`"]

(********************************************************************) 
(********************************************************************) 

(*=========================Simulation==============================*)
ComputeMA=
  Compile[{{theta, _Real, 1},{a,_Real,1}},
    With[{q=Length[theta]},Partition[a,q+1,1].Prepend[-theta,1.0]]];

SimulateAR1[\[Phi]_Real /; Abs[\[Phi]] < 1, n_Integer /; n > 0] := 
  Module[{z1, a = RandomArray[NormalDistribution[0, 1], n - 1]}, 
   z1 = Random[NormalDistribution[0.0, 1.0]]/Sqrt[1.0 - \[Phi]^2]; 
    FoldList[\[Phi]*#1 + #2 & , z1, a]];

SimulateGaussianIMA[\[Theta]_Number, n_Integer /; n > 0, Z1_Number] := 
  Module[{w = SimulateGaussianMA[\[Theta], n - 1]}, FoldList[Plus, Z1, w]];

SimulateGaussianAR1[n_Integer /; n > 0, \[Phi]_ /; Abs[\[Phi]] < 1] := 
  Module[{a = RandomArray[NormalDistribution[0.0, 1.0], n - 1]}, 
   z1 = Random[NormalDistribution[0.0, 1.0]]/Sqrt[1 - \[Phi]^2]; 
    FoldList[\[Phi]*#1 + #2 & , z1, a]];
 
SimulateGaussianMA[(\[Theta]_)?NumberVectorQ, n_Integer /; n > 0] := 
  ComputeMA[N[\[Theta]], RandomArray[NormalDistribution[0.0, 1.0], 
    n + Length[\[Theta]]]];
 
SimulateGaussianMA[\[Theta]_Number, n_Integer /; n > 0] := 
  SimulateGaussianMA[{\[Theta]}, n];

SimulateGaussianAR[(phi_)?InvertibleQ, n_Integer /; n > 0] := 
  Module[{p = Length[phi], a, zp}, 
   a = RandomArray[NormalDistribution[0, 1], n]; If[p == 0, Return[a]]; 
    zp = Random[MultinormalDistribution[Table[0, {p}], 
       Toeplitz[TacvfARMA[phi, {}, p - 1]]]]; 
    Flatten[Prepend[Transpose[FoldList[Prepend[Drop[#1, -1], 
           phi . #1 + #2] & , zp, Drop[a, p - 1]]][[1]], Drop[zp, 1]]]];
 
SimulateGaussianAR[phi_Number, n_Integer /; n > 0] := 
  SimulateGaussianAR[{phi}, n];

SimulateGaussianARMA[(phi_)?InvertibleQ, (theta_)?NumberVectorQ, 
   n_Integer /; n > 0] := 
  Module[{p = Length[phi], q = Length[theta], a, zp, x}, 
   a = RandomArray[NormalDistribution[0, 1], n + q]; 
    If[p == 0 && q == 0, Return[Take[a, n]]]; x = ComputeMA[N[theta], a]; 
    If[p == 0, Return[x]]; zp = 
     Random[MultinormalDistribution[Table[0, {p}], 
       Toeplitz[TacvfARMA[phi, theta, p - 1]]]]; 
    Flatten[Prepend[Transpose[FoldList[Prepend[Drop[#1, -1], 
           phi . #1 + #2] & , zp, Drop[x, p - 1]]][[1]], Drop[zp, 1]]]];
 
SimulateGaussianARMA[phi_Number, (theta_)?NumberVectorQ, 
   n_Integer /; n > 0] := SimulateGaussianARMA[{phi}, theta, n];
 
SimulateGaussianARMA[phi_Number, theta_Number, n_Integer /; n > 0] := 
  SimulateGaussianARMA[{phi}, {theta}, n];
 
SimulateGaussianARMA[(phi_)?InvertibleQ, theta_Number, n_Integer /; n > 0] := 
  SimulateGaussianARMA[phi, {theta}, n];

SimulateARIMA[phi_?NumberVectorQ, theta_?NumberVectorQ, n_Integer] := 
    SimulateGP[TacvfARMA[phi, theta, n - 1]];
 
SimulateARIMA[phi_?NumberVectorQ, theta_?NumberVectorQ, sphi_?NumberVectorQ, 
      stheta_?NumberVectorQ, s_?IntegerQ, n_?IntegerQ, (opts___)?OptionQ] := 
    Module[{fd, diff, gphi, gtheta, nwarm},
      fd = FractionalDifference /. {opts} /. Options[SimulateARIMA]; 
      diff = Differencing /. {opts} /. Options[SimulateARIMA];
      nwarm = WarmUp /. {opts} /. Options[SimulateARIMA]; 
      {gphi, gtheta} = SARMAToARMA[phi, theta, sphi, stheta, s];
      w = SimulateGP[ TacvfARMA[phi, theta, n + nwarm - 1, FractionalDifference -> fd]];
      w = Which[Length[diff] == 0, w, 
          Length[diff] == 2, TSIntegrate[w, diff[[1]], diff[[2]]],
          Length[diff] == 4, 
          TSIntegrate[w, diff[[1]], diff[[2]], diff[[3]], diff[[4]]]];
      Take[w, -n]] ;


SimulateARIMA[(phi_)?NumberVectorQ, (theta_)?NumberVectorQ, (n_)?IntegerQ, 
     (opts___)?OptionQ] := Module[{fd, diff, nwarmup}, 
     fd = FractionalDifference /. {opts} /. Options[SimulateARIMA]; 
      diff = Differencing /. {opts} /. Options[SimulateARIMA]; 
      nwarm = WarmUp /. {opts} /. Options[SimulateARIMA]; 
      SimulateARIMA[phi, theta, {}, {}, 1, n, FractionalDifference -> fd, 
       Differencing -> diff, WarmUp -> nwarm]];

SimulateGP[c_?NumberVectorQ/;Length[c]>1] := 
  Module[{a, z, phik, phi1, vk, phikm1, vkm1, d}, 
   d = Drop[c, 1]; phik = {phi1 = c[[2]]/c[[1]]}; 
    z = RandomArray[NormalDistribution[0, 1], Length[c]]; 
    vk = c[[1]]*(1 - phi1^2); z[[1]] = Sqrt[c[[1]]]*z[[1]]; 
    z[[2]] = Sqrt[vk]*z[[2]] + phi1*z[[1]]; 
    Do[vkm1 = vk; phikm1 = phik; a = DLa[phikm1, d, k - 1]/vk; 
      phik = DLphik[phikm1, a]; vk = vkm1*(1 - a^2); 
      (z[[k]] = DLz[z, phik, vk]; ), {k, 3, Length[c]}]; z];

SimulateTent[n_Integer /; n > 0, z0_:0.177221] := 
  If[z0 < 1 && z0 > 0, N[With[{XSTART = Rationalize[z0]}, 
     NestList[If[#1 < 1/2, 2*#1, 2 - 2*#1] & , XSTART, n]]]];

SimulateIA[(n_)?IntegerQ, (T_)?IntegerQ, (omega_)?VectorQ, (delta_)?VectorQ, 
     (phi_)?VectorQ, (theta_)?VectorQ, (opts___)?OptionQ] := 
    Module[{type}, type = IAType /. {opts} /. Options[SimulateIA]; 
      dif = Differencing /. {opts} /. Options[SimulateIA]; 
      sigsq = SigmaSquared /. {opts} /. Options[SimulateIA]; 
      Sqrt[sigsq]*SimulateARIMA[phi, theta, n, Differencing -> dif] + 
       ICalc[n, T, omega, delta, IAType -> type]];
 
SimulateIA[(n_)?IntegerQ, (T_)?IntegerQ, (omega_)?NumericQ, 
     (delta_)?NumericQ, (phi_)?NumericQ, (theta_)?NumericQ, 
     (opts___)?OptionQ] := Module[{type, dif}, 
     type = IAType /. {opts} /. Options[SimulateIA]; 
      dif = Differencing /. {opts} /. Options[SimulateIA]; 
      sigsq = SigmaSquared /. {opts} /. Options[SimulateIA]; 
      SimulateIA[n, T, {omega}, {delta}, {phi}, {theta}, IAType -> type, 
       Differencing -> dif, SigmaSquared -> sigsq]];
 
SimulateIA[(n_)?IntegerQ, (T_)?VectorQ, (omega_)?ListQ, (delta_)?ListQ, 
     (phi_)?VectorQ, (theta_)?VectorQ, (opts___)?OptionQ] := 
    Module[{type, op, a, z = Array[0 & , n], dif}, 
     If[Length[Union[Length /@ {T, omega, delta}]] > 1, 
       Message[SimulateIA::argserr, Length[T]]]; 
      dif = Differencing /. {opts} /. Options[SimulateIA]; 
      a = IAType /. {opts} /. Options[SimulateIA]; 
      sigsq = SigmaSquared /. {opts} /. Options[SimulateIA]; 
      op = If[Length[a] == Length[T], a, Flatten[Table[a, {Length[T]}]]]; 
      Do[z += Sqrt[sigsq]*SimulateARIMA[phi, theta, n, Differencing -> dif] + 
         ICalc[n, T[[i]], omega[[i]], delta[[i]], IAType -> op[[i]]], 
       {i, Length[T]}]; z];

(*=========================SDF==============================*)
SdfPlot[(z_)?NumberVectorQ, (q_)?IntegerNumberVectorQ, (opts___)?OptionQ] := 
  Module[{plotit,norq,pr,pl,ylog, xlog, p, f, ms, ic, alg, res, ar, x, n = Length[z]}, 
   xlog = LoggedF /. {opts} /. Options[SdfPlot]; 
    ylog = LoggedSdf /. {opts} /. Options[SdfPlot]; 
    norq = Normalized /. {opts} /. Options[SdfPlot]; 
    plotit = Plot /. {opts} /. Options[SdfPlot]; 
    pr = PlotRange /. {opts} /. Options[SdfPlot]; 
    pl = PlotLabel /. {opts} /. Options[SdfPlot];
    delta = SamplingFrequency/.{opts}/.Options[SdfPlot]; 
    optListPlot = FilterOptions[ListPlot, opts];
    {f, p} = If[ListQ[q], 
        {f, p} = If[norq,NormalizedPeriodogram[z],Periodogram[z]]; 
       p = If[Length[q] == 0, p, Fold[DaniellSmoother[#1, #2] & , p, q]]; 
       {f, p}, ms = ModelSelection /. {opts} /. Options[SdfPlot]; 
       alg = ARAlgorithm /. {opts} /. Options[SdfPlot]; 
       ic = Which[MatchQ[ms, False], False, MatchQ[ms, None], False, 
         MatchQ[ms, Automatic], True, MatchQ[ms, "AIC"], "AIC", 
         MatchQ[ms, "BIC"], "BIC", True, True]; 
       res = If[MatchQ[alg, "YuleWalker"], 
         PacfDL[Sacvf[z, MaxLag -> q], LinearPredictor -> True], 
         PacfBurg[z, q, LinearPredictor -> True]]; 
       qopt = If[StringQ[ic], Which[MatchQ[ic, "BIC"], 
          ICFromPacf[n, Pacf /. res, Criterion -> "BIC"]], 
         If[ic, ICFromPacf[n, Pacf /. res, Criterion -> "AIC"], q]]; 
       res = If[MatchQ[alg, "YuleWalker"], 
         PacfDL[Sacvf[z, MaxLag -> qopt], LinearPredictor -> True], 
         PacfBurg[z, qopt, LinearPredictor -> True]]; 
       ar = ARCoefficients /. res; resvar = ResidualVariance /. res; 
       x = SdfARMA[ar, {}, BenchmarkTimeSeriesLength -> Length[z],
           Normalized->norq]; 
       MapAt[resvar*#1 & , Transpose[x], 2]]; f = If[xlog, Log[2, f], f]; 
    p = If[ylog, Log[2, p], p]; 
    If[plotit,
     ListPlot[Transpose[{f*delta, p}], Axes -> False, Frame -> True, 
        PlotRange -> pr, PlotJoined -> True, PlotLabel->pl, optListPlot],
     Transpose[{f*delta, p}]
    ] ];
 
NPSdfPlot[(z_)?NumberVectorQ,
    (opts___)?OptionQ] := 
  Module[{plotit,norq,pr,pl,ylog, xlog, p, f, q}, 
    q = SdfSpan /. {opts} /. Options[NPSdfPlot]; 
    xlog = LoggedF /. {opts} /. Options[NPSdfPlot]; 
    ylog = LoggedSdf /. {opts} /. Options[NPSdfPlot]; 
    norq = Normalized /. {opts} /. Options[NPSdfPlot]; 
    plotit = Plot /. {opts} /. Options[NPSdfPlot]; 
    pr = PlotRange /. {opts} /. Options[NPSdfPlot]; 
    pl = PlotLabel /. {opts} /. Options[NPSdfPlot]; 
    {f, p} = If[norq,NormalizedPeriodogram[z],Periodogram[z]]; 
    p = Which[
          q==Automatic, Print["ok use Lee"];LeeSmoother[p],
          Length[q] == 0, p, 
          True, Print["ok-Daniell"];Fold[DaniellSmoother[#1, #2] & , p, q]
        ]; 
    f = If[xlog, Log[2, f], f]; 
    p = If[ylog, Log[2, p], p]; 
    If[plotit,
     ListPlot[Transpose[{f, p}], Axes -> False, Frame -> True, 
        PlotRange -> pr, PlotJoined -> True, PlotLabel->pl],
     Transpose[{f, p}]
    ] ];
 
ARSdfPlot[(z_)?NumberVectorQ, (q_)?IntegerNumberVectorQ, (opts___)?OptionQ] := 
  Module[{plotit,norq,pr,pl,ylog, xlog, p, f, ms, ic, alg, 
          res, ar, x, n = Length[z], delta}, 
    xlog = LoggedF /. {opts} /. Options[ARSdfPlot]; 
    ylog = LoggedSdf /. {opts} /. Options[ARSdfPlot]; 
    norq = Normalized /. {opts} /. Options[ARSdfPlot]; 
    plotit = Plot /. {opts} /. Options[ARSdfPlot]; 
    pr = PlotRange /. {opts} /. Options[ARSdfPlot]; 
    pl = PlotLabel /. {opts} /. Options[ARSdfPlot]; 
    ms = ModelSelection /. {opts} /. Options[ARSdfPlot]; 
    delta = SamplingFrequency /. {opts} /. Options[ARSdfPlot]; 
    alg = EstimationAlgorithm /. {opts} /. Options[ARSdfPlot];
    optListPlot = FilterOptions[ListPlot, opts]; 
    ic = Which[MatchQ[ms, "AIC"], 1, 
               MatchQ[ms, "BIC"], 2, 
               True, 3]; 
    res = If[MatchQ[alg, "YuleWalker"], 
         PacfDL[Sacvf[z, MaxLag -> q], LinearPredictor -> True], 
         PacfBurg[z, q, LinearPredictor -> True]]; 
    qopt = Which[ic==1, ICFromPacf[n, Pacf /. res, Criterion -> "AIC"], 
                 ic==2, ICFromPacf[n, Pacf /. res, Criterion -> "BIC"],
                 True,  q];
    res = If[ms===Automatic,PacfBayes[z,q,LinearPredictor->True],
            If[MatchQ[alg, "YuleWalker"], 
             PacfDL[Sacvf[z, MaxLag -> qopt], LinearPredictor -> True], 
             PacfBurg[z, qopt, LinearPredictor -> True]]]; 
    ar = ARCoefficients /. res; resvar = ResidualVariance /. res; 
    x = SdfARMA[ar, {}, BenchmarkTimeSeriesLength -> Length[z],
           Normalized->norq]; 
    {f, p} = If[norq, Transpose[x], MapAt[resvar*#1 & , Transpose[x], 2]]; 
    f = delta*f;
    f = If[xlog, Log[2, f], f]; 
    p = If[ylog, Log[2, p], p];
    If[plotit,
     ListPlot[Transpose[{f, p}], Axes -> False, Frame -> True, 
        PlotRange -> pr, PlotJoined -> True, PlotLabel->pl,optListPlot],
     Transpose[{f, p}]
    ] ];
 
NormalizedPeriodogram[z_?NumberVectorQ]:=With[{v=VarianceMLE[z]},
MapAt[Divide[#,v]&,Periodogram[z],2]];

Periodogram=
  Compile[{{z,_Real,1}},
    With[{f=Range[Floor[Length[z]/2]]/Length[z]},{f,
        Take[Drop[With[{t=Fourier[z]},t*Conjugate[t]//Re],1],Length[f]]}]];

DaniellSmoother[(p_)?VectorQ, q_Integer] := 
  Partition[Join[Reverse[Take[Drop[p, 1], q]], p, 
      Reverse[Take[Drop[p, -1], -q]]], 2*q + 1, 1] . 
    Append[Prepend[Table[1/(2*q), {2*q - 1}], 1/(4*q)], 1/(4*q)] /; q > 0;
 
CumulativePeriodogramTest[z_] := 
  Module[{f, p, cp, K = 1.629999999999999/Sqrt[Floor[(Length[z] - 1)/2]]}, 
   {f, p} = Periodogram[z]; cp = Drop[FoldList[Plus, 0, p], 1]/Plus @@ p; 
    ListPlot[Transpose[{f, cp}], Axes -> False, Frame -> True, 
     Epilog -> {{RGBColor[0, 1, 0], Line[{{0, 0}, {0.5, 1}}]}, 
       {RGBColor[1, 0, 0], Line[{{K/2, 0}, {0.5, 1 - K}}], 
        Line[{{0, K}, {0.5, 1 + K}}]}}]];
 
IteratedDaniellWeights[(q_)?VectorQ] := 
  Module[{pulse = Table[0, {2*Plus @@ q + 3}]}, 
   pulse[[1 + (Length[pulse] - 1)/2]] = 1; 
    Fold[DaniellSmoother[#1, #2] & , pulse, q]];
 
OptimalPeriodogramSmoother[z_?NumberVectorQ] := 
  Module[{rulesQ,n = Length[z], f, p, per, sumsqI}, 
   per = {f, p} = Periodogram[z]; 
   sumsqI = 2*Plus @@ (p^2); 
   EstimatedRisk[q_] := 
                2*Plus @@ ((DaniellSmoother[p, q] - p)^2)/n - 
                   (1 - 2/(2*q))/(2*n)*sumsqI;
   IndexOfMinimum[Table[EstimatedRisk[q], 
                         {q, 1, Ceiling[10*Log[10, n]]}]]
      ]; 

CISdf[z_?NumberVectorQ, (opts___)?OptionQ] := 
CISdf[z,{OptimalPeriodogramSmoother[z]},opts];

CISdf[z_?NumberVectorQ, q_List, (opts___)?OptionQ] := 
  Module[{f, p, chisq, a, norq}, 
    norq = Normalized /. {opts} /. Options[CISdf];
    {f, p} = If[norq,NormalizedPeriodogram[z],Periodogram[z]]; 
    p = If[Length[q] == 0, p, Fold[DaniellSmoother[#1, #2] & , p, q]]; 
    df = 2./Plus @@ (IteratedDaniellWeights[q]^2); 
    chisq = Quantile[ChiSquareDistribution[df], {0.975, 0.025}]; 
    a = df/chisq; Prepend[Transpose[Outer[Divide, p, a]], f]];
 
SdfCheck[(z_)?NumberVectorQ, (sdf_)?MatrixQ, (opts___)?OptionQ] := 
SdfCheck[z, sdf,{OptimalPeriodogramSmoother[z]}, opts];

SdfCheck[(z_)?NumberVectorQ, (sdf_)?MatrixQ, q_List, (opts___)?OptionQ] := 
  Module[{pl,xlog, ylog, xf, yf, sdf2, c, smax, smin}, 
    xlog = LoggedF /. {opts} /. Options[SdfCheck]; 
    ylog = LoggedSdf /. {opts} /. Options[SdfCheck]; 
    pl = PlotLabel /. {opts} /. Options[SdfCheck]; 
    xf = Which[xlog, Log[2, #1] & , True, #1 & ]; 
    yf = Which[ylog, Log[2, #1] & , True, #1 & ]; 
    sdf2 = Transpose[MapAt[yf, Transpose[sdf], {2}]]; 
    sdf2 = Transpose[MapAt[xf, Transpose[sdf2], {1}]]; 
    c = CISdf[z, q, Normalized->True]; 
    c = MapAt[yf, c, {{2}, {3}}]; c = MapAt[xf, c, {1}]; 
    csdf = {c, Transpose[sdf2]}; smax = Max[(Max /@ Drop[#1, 1] & ) /@ csdf]; 
    smin = Min[(Min /@ Drop[#1, 1] & ) /@ csdf]; 
    ListPlot[sdf2, PlotLabel->pl,PlotRange -> {Automatic, {smin, smax}}, Axes -> False, 
     PlotJoined -> True, Frame ->True,  
     Epilog -> {RGBColor[1, 0, 1], Line[Transpose[{First[c], c[[2]]}]], 
       Line[Transpose[{First[c], c[[3]]}]]}]];

SdfARMAPlot[(phi_)?NumberVectorQ, (theta_)?NumberVectorQ, 
   (opts___)?OptionQ] := 
  Module[{pr, pl, sdf, framelab, xlog, ylog, d, f, n}, 
   xlog = LoggedF /. {opts} /. Options[SdfARMAPlot]; 
    ylog = LoggedSdf /. {opts} /. Options[SdfARMAPlot]; 
    d = FractionalDifference /. {opts} /. Options[SdfARMAPlot]; 
    n = BenchmarkTimeSeriesLength /. {opts} /. Options[SdfARMAPlot]; 
    pr = PlotRange /. {opts} /. Options[SdfARMAPlot]; 
    pl = PlotLabel /. {opts} /. Options[SdfARMAPlot]; 
    {f, sdf} = Transpose[SdfARMA[phi, theta, FractionalDifference -> d, 
       TimeSeriesLength -> n]]; f = If[xlog, Log[2, f], f]; 
    sdf = If[ylog, Log[2, sdf], sdf]; 
    framelab = Which[ !(xlog || ylog), {"frequency", "sdf"},  !xlog && ylog, 
      {"frequency", "Log[2,sdf]"}, xlog && ylog, 
      {"Log[2,frequency]", "Log[2,sdf]"}, True, {"Log[2,frequency]", "sdf"}]; 
    ListPlot[Transpose[{f, sdf}], PlotRange -> pr, PlotLabel->pl, 
     PlotJoined -> True, 
     Frame -> True, Axes -> False, FrameLabel -> framelab, 
     RotateLabel -> False]];

SdfARMA[(phi_)?NumberVectorQ, (theta_)?NumberVectorQ, (opts___)?OptionQ] := 
  Module[{d, sdf, v, norq}, 
    d = FractionalDifference /. {opts} /. Options[SdfARMA]; 
    v = InnovationVariance /. {opts} /. Options[SdfARMA]; 
    n = BenchmarkTimeSeriesLength /. {opts} /. Options[SdfARMA]; 
    norq = Normalized /. {opts} /. Options[SdfARMA]; 
    sdf = Drop[Which[Length[phi] == 0 && Length[theta] == 0, 
       Table[1, {512}], Length[phi] == 0, SdfMA[theta], Length[theta] == 0, 
       SdfAR[phi], True, SdfAR[phi]*SdfMA[theta]], Ceiling[1024/n]];
    v=If[norq,1.0/First[TacvfARMA[phi,theta,1,FractionalDifference->d]],v];
    Transpose[{Range[Ceiling[1024/n], 511]/1024., 
       v*If[d == 0, sdf, sdf*SdfFDWN[d, n]]}]];
 
SdfMA[(theta_)?NumberVectorQ] := 
  With[{pFFT = 9}, (2^(pFFT + 1)*Re[#1*Conjugate[#1]] & )[Take[
     Fourier[Join[{1}, -theta, 
       Table[0, {2^(pFFT + 1) - 1 - Length[theta]}]]], 2^pFFT]]];
 
SdfAR[(phi_)?NumberVectorQ] := 
  With[{pFFT = 9}, 1./
    (2^(pFFT + 1)*Re[#1*Conjugate[#1]] & )[Take[Fourier[Join[{1}, -phi, 
        Table[0, {2^(pFFT + 1) - 1 - Length[phi]}]]], 2^pFFT]]];

SdfFDWN=Compile[{{d, _Real},{n, _Integer}},
    With[{m=512},
    ((2.0*Sin[Pi*#])^(-2*d))&/@
        (Range[Ceiling[(2*m)/n],m-1]/(2*m)//N)]];

SdfARMAvsFARMAPlot[(armaphi_)?NumberVectorQ, (armatheta_)?NumberVectorQ, 
   (farmaphi_)?NumberVectorQ, (farmatheta_)?NumberVectorQ, 
   d_Real /; Abs[d] < 0.5, (opts___)?OptionQ] := 
  Module[{sdfarma, sdffarma, framelab, xlog, ylog, 
          f, n,pl,pr,plegend,pstyle}, 
   xlog = LoggedF /. {opts} /. Options[SdfARMAvsFARMAPlot]; 
    ylog = LoggedSdf /. {opts} /. Options[SdfARMAvsFARMAPlot]; 
    n = BenchmarkTimeSeriesLength /. {opts} /. Options[SdfARMAvsFARMAPlot]; 
    pr = PlotRange /. {opts} /. Options[SdfARMAvsFARMAPlot]; 
    pl = PlotLabel /. {opts} /. Options[SdfARMAvsFARMAPlot]; 
    plegend = PLegend /. {opts} /. Options[SdfARMAvsFARMAPlot]; 
    pstyle = PStyle /. {opts} /. Options[SdfARMAvsFARMAPlot]; 
    {f, sdfarma} = 
     Transpose[SdfARMA[armaphi, armatheta, TimeSeriesLength -> n, 
               Normalized->True]]; 
    {f, sdffarma} = 
     Transpose[SdfARMA[farmaphi, farmatheta, TimeSeriesLength -> n, 
       FractionalDifference -> d,Normalized->True]]; 
    f = If[xlog, Log[2, f], f]; 
    {sdfarma, sdffarma} = 
     If[ylog, (Log[2, #1] & ) /@ {sdfarma, sdffarma}, {sdfarma, sdffarma}]; 
    framelab = Which[ !(xlog || ylog), {"frequency", "sdf"},  !xlog && ylog, 
      {"frequency", "Log[2,sdf]"}, xlog && ylog, 
      {"Log[2,frequency]", "Log[2,sdf]"}, True, {"Log[2,frequency]", "sdf"}]; 
    MultipleListPlot[Transpose[{f, sdfarma}], Transpose[{f, sdffarma}], 
     PlotRange -> pr, PlotJoined -> {True, True}, Frame -> True, 
     Axes -> False, FrameLabel -> framelab, RotateLabel -> False, 
     PlotStyle -> pstyle,
     SymbolShape -> {None, None}, PlotLabel->pl,
     PlotLegend -> plegend]];
 
(*===========================Visualization==============================*)
AcfPlot[c_, (opts___)?OptionQ] :=
   Module[{acfrange, acfticks,c2=Prepend[c,1]},
   acfrange=If[MemberQ[Map[(#>0)&, c2], False],{-1, 1}, {0, 1}];
   frameticks={Range[0,Length[c2],Max[1,Round[Length[c2]/10.0]]], 
             Range[acfrange[[1]],acfrange[[2]],.5]}//Chop;
   Show[Graphics[
    {RGBColor[0,0,1], Thickness[0.01],
     Map[Line, 
         Map[{{#1 - 1, 0}, {#1 - 1, c2[[#1]]}} & , Range[Length[c2]]]
        ]
    },
                opts,
                Frame->True,
                Axes->False,
                AxesOrigin->{0,0}, 
                FrameLabel->{"Lag", "Acf"},
                PlotRange->acfrange,
                FrameTicks->frameticks,
                GridLines->{{}, Range[-1,1,.5]}
       ]]
	];

NPlot[(data_)?NumberVectorQ, (plotopts___)?OptionQ] := 
  ListPlot[Transpose[{Table[Quantile[NormalDistribution[0, 1], 
       (i - 0.5)/Length[data]], {i, Length[data]}], Sort[data]}], 
   PlotRange -> All, Prolog -> {PointSize[0.02]}, 
   Epilog -> quartileLine[data], plotopts];
 
quartileLine[(x_)?NumberVectorQ] := 
  Module[{slope, intercept, quartilesX, quartilesY, extremesX}, 
   quartilesY = {Quantile[x, 0.25], Quantile[x, 0.75]}; 
    quartilesX = 
     {Quantile[NormalDistribution[0, 1], 0.25], 
      Quantile[NormalDistribution[0, 1], 0.75]}; 
    slope = (quartilesY[[2]] - quartilesY[[1]])/
      (quartilesX[[2]] - quartilesX[[1]]); 
    intercept = quartilesY[[1]] - slope*quartilesX[[1]]; 
    extremesX = 
     {Quantile[NormalDistribution[0, 1], 0.5/Length[x]], 
      Quantile[NormalDistribution[0, 1], 1 - 0.5/Length[x]]}; 
    Line[{{extremesX[[1]], intercept + slope*extremesX[[1]]}, 
      {extremesX[[2]], intercept + slope*extremesX[[2]]}}]];
 
 
BWPlot[(data_)?NumberVectorMatrixQ, (plotopts___)?OptionQ] := 
  Module[{datadim, k, datapts, whiskers, box, outsidepts, medpt, dmax, dmin, 
    coldata, Q1, Q3, uplim, dnlim, outside, jitter, drange, epsilon, 
    boxwidth = 0.4}, datadim = Dimensions[data]; 
    k = If[Length[datadim] == 2, datadim[[2]], 1]; 
    datapts = outsidepts = whiskers = box = medpt = dmax = dmin = {}; 
    Do[coldata = If[k == 1, data, Column[data, i]]; 
      datapts = 
       Join[datapts, Transpose[{coldata, Table[i, {Length[coldata]}]}]]; 
      medpt = Join[medpt, {PointSize[0.04], 
         Point[{Quantile[coldata, 0.5], i}]}]; 
      Q1 = Quantile[coldata, 0.25]; Q3 = Quantile[coldata, 0.75]; 
      box = Join[box, {RGBColor[0.690207, 0.7685929999999999, 0.870602], 
         Polygon[{{Q1, i - boxwidth}, {Q1, i + boxwidth}, 
           {Q3, i + boxwidth}, {Q3, i - boxwidth}}]}]; 
      step = 1.5*(Q3 - Q1); uplim = Q3 + step; dnlim = Q1 - step; 
      upadj = Max[Select[coldata, #1 <= uplim & ]]; 
      dnadj = Min[Select[coldata, #1 >= dnlim & ]]; 
      whiskers = 
       Join[whiskers, {Thickness[0.005], Line[{{dnadj, i}, {Q1, i}}], 
         Line[{{Q3, i}, {upadj, i}}]}]; 
      outside = 
       Join[Select[coldata, #1 > uplim & ], Select[coldata, #1 < dnlim & ]];
      jitter = Table[i + (Random[] - 1/2)/6, {Length[outside]}]; 
      outsidepts = 
       Join[outsidepts, (Circle[#1, Offset[{5, 5}]] & ) /@ 
         Transpose[{outside, jitter}]]; dmax = Max[dmax, Max[coldata]]; 
      (dmin = Min[dmin, Min[coldata]]; ), {i, k}]; 
    epsilon = (dmax - dmin)*0.03; drange = {dmin - epsilon, dmax + epsilon}; 
    ListPlot[datapts, Ticks -> {Automatic, None}, 
     PlotStyle -> {AbsolutePointSize[0]}, PlotRange -> {drange, {0, k + 1}}, 
     Axes -> {True, False}, Epilog -> {outsidepts, box, medpt, whiskers}, 
     plotopts]];

BNPlot[(x_)?NumberVectorQ] := 
  Show[GraphicsArray[{{NPlot[x, DisplayFunction -> Identity]}, 
     {BWPlot[x, DisplayFunction -> Identity]}}]];

TimeSeriesPlot[z_?NumberVectorQ, sublength_:200, aspectratio_:0.25, opts_?OptionQ] := 
  If[(sublength>0)&&(aspectratio>0),
   Module[{z2, numpanels, itake, g, eps, optListPlot},
   eps = 0.05*(Max[z]-Min[z]); 
   optListPlot=FilterOptions[ListPlot, opts];
   If[sublength >= Length[z], numpanels = 1, 
     z2 = Partition[Join[Transpose[{Range[Length[z]], z}], 
         Table[0, {sublength}]], sublength, sublength - 1]; 
      numpanels = Length[z2]];  
    If[numpanels > 1, itake = First[First[-1 + Position[Last[z2], 0]]]; 
      z2[[numpanels]] = 
       Join[Drop[z2[[numpanels - 1]], itake], Take[Last[z2], itake]]; 
      g = (ListPlot[#1, PlotJoined -> True, 
           PlotRange -> {All, {Min[z]-eps, Max[z]+eps}}, Frame -> True, 
           Axes -> False, AspectRatio -> 0.25, 
           DisplayFunction -> Identity] & ) /@ z2; 
      Show[GraphicsArray[Partition[g, 1]]], 
     ListPlot[z, PlotJoined -> True, PlotRange -> {All, {Min[z]-eps, Max[z]+eps}}, 
      Frame -> True, Axes -> False, AspectRatio -> 0.25, optListPlot]]]
   ];

AnimateTimeSeriesPlot[(z_)?NumberVectorQ, 
          sublength_Integer /; sublength > 0, 
          advance_Integer/;advance>0, 
          aspectratio_:0.25] := 
    Module[{numpanels, itake, g}, 
     z2 = Partition[Join[Transpose[{Range[Length[z]], z}]], sublength, 
        Min[advance, sublength - 1]]; 
      z2 = Append[z2, Take[Join[Transpose[{Range[Length[z]], z}]], 
         -sublength]]; numpanels = Length[z2]; 
      g = (ListPlot[#1, PlotJoined -> True, PlotRange -> 
           {All, {Min[z], Max[z]}}, Frame -> True, FrameTicks -> 
           {{Min[First[Transpose[#1]]]}, Automatic, None, None}, 
          Axes -> False, AspectRatio -> 0.25, DisplayFunction -> 
           Identity] & ) /@ z2; ShowAnimation[g]];
 
PhasePlot[(z_)?NumberVectorQ, lag_Integer /; lag > 0] := 
  ListPlot[Reverse /@ 
    Drop[Transpose[Reverse[NestList[RotateLeft[#1, lag] & , z, 1]]], -lag], 
   PlotJoined -> True];
 
PlotARMAFactors[phi_] := PlotARMAFactors[phi, {}];

PlotARMAFactors[{phi_List, theta_List}] := PlotARMAFactors[phi, theta];

PlotARMAFactors[phi_, theta_] := 
Module[{framelab, framelabar},
   framelabar = StringForm["phi = ``", phi];
   framelab = FontForm[If[Length[theta]==0, 
                          framelabar, 
                          {framelabar,StringForm["theta = ``", theta]}
                      ],
              {"Helvetica-Bold", 8}];
  Show[
     Graphics[{ 
                 {RGBColor[1,1,0], Disk[{0, 0}, 1, {-Pi, Pi}]},
                 {RGBColor[0,0,1], PointSize[0.02],
                  Map[Point,
                   (Transpose[{Re[#1],Im[#1]}])&[1/ARRoots[phi]]
                     ]
                 },
                 {RGBColor[0.4,0.4,0], PointSize[0.02],
                  Map[Point,
                   (Transpose[{Re[#1],Im[#1]}])&[1/ARRoots[theta]]
                     ]
                 }
               },
                AspectRatio->Automatic,
                Frame->True,
                Axes->True,
                AxesOrigin->{0,0}, 
                AxesLabel->{"Real", "Imaginary"},
                GridLines->{Range[-1,1,.1], Range[-1,1,.1]},
                FrameLabel->{framelab, ""}
              ] 
      ]
];

PlotSARMAFactors[{phi_List, theta_List, PHI_List, THETA_List, S_}] := 
                PlotSARMAFactors[phi, theta, PHI, THETA, S];

PlotSARMAFactors[phi_, theta_, PHI_, THETA_, S_] := 
Module[{framelab},
   framelab = FontForm[ToString[NumberForm[#,2]&/@Flatten[{phi,theta,PHI,THETA,S}]],
              {"Helvetica-Bold", 10}];
  Show[
      Graphics[{ 
                 {RGBColor[1,1,0], Disk[{0, 0}, 1, {-Pi, Pi}]},
                 {RGBColor[0,0,1], PointSize[0.02],
                  Map[Point,
                   (Transpose[{Re[#1],Im[#1]}])&[1/ARRoots[phi]]
                     ]
                 },
                 {RGBColor[0,1,0], PointSize[0.015],
                  Map[Point,
                   (Transpose[{Re[#1],Im[#1]}])&[1/SARRoots[PHI, S]]
                     ]
                 },
                 {RGBColor[0.4,0,0.4], PointSize[0.02],
                  Map[Point,
                   (Transpose[{Re[#1],Im[#1]}])&[1/ARRoots[theta]]
                     ]
                 },
                 {RGBColor[1,0,0], PointSize[0.015],
                  Map[Point,
                   (Transpose[{Re[#1],Im[#1]}])&[1/SARRoots[THETA,S]]
                     ]
                 }
               },
                AspectRatio->Automatic,
                Frame->True,
                Axes->True,
                AxesOrigin->{0,0}, 
                AxesLabel->{"Real", "Imaginary"},
                GridLines->{Range[-1,1,.1], Range[-1,1,.1]},
                FrameLabel->{framelab, ""}
              ] 
      ]
];

(*=========================TACF, etc.==============================*)
getc=Compile[{{p,_Integer},{q,_Integer},{phi,_Real,1},{theta,_Real,1}},
      Which[q\[Equal]0,{1.0},p\[Equal]0,Prepend[-theta,1.0],True,
        First/@FoldList[Prepend[Take[#1,p-1],-#2+phi.#1]&,
            Prepend[Table[0.0,{p-1}],1],theta]]];

getb=Compile[{{theta, _Real, 1},{c, _Real, 1},{s,_Integer},{rp1, _Integer}},
    b=Table[Take[c,s+1-i].Drop[Prepend[theta,-1.0],i-1],{i,s}];
    If[rp1>s,Join[b,Table[0.0,{rp1-s}]],b]];

getg=Compile[{{phi,_Real, 1},{b,_Real,1},{p, _Integer},{rp1, _Integer}},
    If[p==0,-b,
      LinearSolve[
        Table[With[{ij=i-j},
              Which[ij==0,-1,(ij<0)||(ij>p),0,True,
                phi\[LeftDoubleBracket]ij\[RightDoubleBracket]]]+
            If[j==1,0.0,
              With[{ij=i+j-2},
                Which[ij==0,-1.0,(ij<0)||(ij>p),0,True,
                  phi\[LeftDoubleBracket]ij\[RightDoubleBracket]]]],{i,rp1},{
            j,rp1}],b]]];

getr=Compile[{{phi,_Real,1},{g,_Real,1},{p,_Integer},{rp1,_Integer},{
        MaxLag,_Integer}},
    If[MaxLag<rp1,Take[g,MaxLag+1],
      If[p==0,Join[g,Table[0.0,{MaxLag-rp1+1}]],
        Join[Drop[g,-1],
          First/@NestList[Take[Prepend[#,phi.#],p]&,Take[Reverse[g],p],
              MaxLag-rp1+1]]]]];

TacvfARMA[(phi_List)?NumberVectorQ,(theta_List)?NumberVectorQ,
      MaxLag_Integer,(opts___)?OptionQ]:=
 Module[{q,p,A,b,c,phifun,rp1,qp1, lagtrunc,lagtrunc2,d,Gphi,Gtheta2},
  If[And @@ InvertibleQ /@ {phi, {}}//Not, Message[TacvfARMA::singular];Return[{}]];
  lagtrunc2=LagTruncation/.{opts}/.Options[TacvfARMA]; 
  d=FractionalDifference/.{opts}/.Options[TacvfARMA]; 
(*next 2 lines added to adhere to new requirement with Compile function *)
(*but there is obviously a better way to do this! *)
  Gphi=If[Length[phi]==0,{0.0},phi];
  Gtheta=If[Length[theta]==0,{0.0},theta];
  p=Length[Gphi];
  q=Length[Gtheta];
  rp1=1+Max[p,q]; 
  qp1=1+q; 
  Which[q\[Equal]0&&p\[Equal]0&&d\[Equal]0,
        Prepend[Table[0.0,{MaxLag}],1],
        q\[Equal]0&&p\[Equal]0&&d\[NotEqual]0,
        TacvfFDWN[d,MaxLag],
        True,
        c=getc[p,q,Gphi,Gtheta]; 
          b=getb[Gtheta,c,qp1,rp1]; 
          g=getg[Gphi,b,p,rp1]; 
          Which[d\[Equal]0,
                getr[Gphi,g,p,rp1,MaxLag],
                True,
                With[{lagtrunc=Max[lagtrunc2,2^Ceiling[Log[2.,MaxLag]]]},
                    Take[Reverse[Chop[Sqrt[
                     2*lagtrunc+1]*Drop[InverseFourier[MapThread[Times,(Fourier[#1]&)/@(Join[Reverse[#1],Drop[#1,1]]&)/@
                     {getr[Gphi,g,p,rp1, lagtrunc],TacvfFDWN[d,lagtrunc]}]],
                      lagtrunc]]],MaxLag+1]]
         ]
  ]
];
 
TacvfARMA[phi_List, theta_List, MaxLag_Integer, (opts___)?OptionQ] := 
    Module[{q = Length[theta], p = Length[phi], A, b, c, theta2, phifun, rp1, 
          qp1, d}, 
        d = FractionalDifference /. {opts} /. Options[TacvfARMA]; 
        If[d != 0, Message[TacvfARMA::"unavail"], rp1 = 1 + Max[p, q]; 
          qp1 = 1 + q; 
          If[q == 0 && p == 0, Return[Prepend[Table[0, {MaxLag}], 1]]]; 
          c = Which[q == 0, {1}, p == 0, Prepend[-theta, 1], True, 
              First /@ 
                FoldList[Prepend[Take[#1, p - 1], -#2 + phi.#1] &, 
                  Prepend[Table[0, {p - 1}], 1], theta]]; 
          theta2 = Prepend[theta, -1]; 
          b = Table[Take[c, qp1 + 1 - i].Drop[theta2, i - 1], {i, qp1}]; 
          b = If[rp1 > qp1, Join[b, Table[0, {rp1 - qp1}]], b]; 
          g = If[p == 0, -b, 
              A = Table[
                  With[{ij = i - j}, 
                      Which[ij == 0, -1, ij < 0 || ij > p, 0, True, 
                        phi[[ij]]]] + 
                    If[j == 1, 0, 
                      With[{ij = i + j - 2}, 
                        Which[ij == 0, -1, ij < 0 || ij > p, 0, True, 
                          phi[[ij]]]]], {i, rp1}, {j, rp1}]; 
              Check[LinearSolve[A, b], {}]]; 
          If[SameQ[g, {}], Message[TacvfARMA::"singular"]; Return[g],
            
            If[MaxLag < rp1 - 1, Take[g, MaxLag + 1], 
              If[p == 0, Join[g, Table[0, {MaxLag - rp1 + 1}]], 
                Join[Drop[g, -1], 
                  First /@ 
                    NestList[Take[Prepend[#1, phi.#1], p] &, 
                      Take[Reverse[g], p], MaxLag - rp1 + 1]]]]]]] /; 
      MaxLag >= 0;

CovarianceMatrixARMA[phi_List, theta_List, n_Integer, opts___] := 
    Toeplitz[TacvfARMA[phi, theta, n-1, opts]] /;n>0 ;

InverseCovarianceMatrixARMA[phi_List, theta_List, n_Integer, 
     (opts___)?OptionQ] := 
 Module[
  {theta2, i, j, k, WMInv, E1, E2, W, A1, A2, A3, A22, B, DInv, Omega1, 
   Omega2, dt, S1, S2, S3, S22, T1, T2, T3, T22, icov, ifdet,
   p = Length[phi], q = Length[theta], 
   \[GothicCapitalD], \[GothicCapitalC], \[DoubleStruckL] },
   ifdet = Determinant /. {opts} /. Options[InverseCovarianceMatrixARMA]; 
   Which[2*Max[p,q] >= n,
           With[{icov=CovarianceMatrixARMA[phi,theta,n]},
             If[ifdet,{Det[icov],Inverse[icov]},Inverse[icov]]],
         p == 0 && q == 0, 
           If[ifdet,{1,IdentityMatrix[n]}, IdentityMatrix[n]], 
         q == 0, 
           InverseCovarianceMatrixAR[phi, n, Determinant -> ifdet], 
         p == 0, 
           InverseCovarianceMatrixMA[theta, n, Determinant -> ifdet], 
         True, 
           \[DoubleStruckL] = Max[p, q]; 
           theta2 = Join[{-1}, theta, Table[0, {\[DoubleStruckL] - q}]]; 
           B = Table[0, {\[DoubleStruckL]}, {\[DoubleStruckL]}]; 
           For[j = 1, j <= \[DoubleStruckL], ++j, 
             For[i = 1, i <= \[DoubleStruckL], ++i, 
               B[[i,j]] = Which[
                 j > i, B[[j,i]], 
                 j <= i, 
                   Sum[theta2[[1 + \[DoubleStruckL] + k - i]]*
                     theta2[[1 + \[DoubleStruckL] + k - j]], 
                     {k, 1, Min[i, j]}]
    ]]]; 
    B = InvertRC[B]; WMInv = SubMatrix[CovarianceMatrixARMA[{}, theta, 
             n], {1, 1}, {\[DoubleStruckL], \[DoubleStruckL]}] - B; 
    E1 = SubMatrix[InverseCovarianceMatrixAR[phi, n, Determinant -> 
               False] . CovarianceMatrixARMA[phi, theta, n], {1, 1}, 
            {\[DoubleStruckL], \[DoubleStruckL]}] - WMInv; 
    E2 = SubMatrix[InverseCovarianceMatrixAR[theta, n, Determinant -> 
               False] . CovarianceMatrixARMA[theta, phi, n] - 
            InverseCovarianceMatrixAR[phi, n, Determinant -> False], {1, 1}, 
           {\[DoubleStruckL], \[DoubleStruckL]}]; 
    W = CovarianceMatrixARMA[theta, {}, n]; 
    A1 = SubMatrix[W, {1, 1}, {\[DoubleStruckL], \[DoubleStruckL]}]; 
    A2 = SubMatrix[W, {\[DoubleStruckL] + 1, 1}, 
           {n - 2*\[DoubleStruckL], \[DoubleStruckL]}]; 
    A3 = SubMatrix[W, {n - \[DoubleStruckL] + 1, 1}, {\[DoubleStruckL], 
            \[DoubleStruckL]}]; A22 = SubMatrix[W, {\[DoubleStruckL] + 1, 
            \[DoubleStruckL] + 1}, {n - 2*\[DoubleStruckL], 
            n - 2*\[DoubleStruckL]}]; \[GothicCapitalD] = 
          IdentityMatrix[\[DoubleStruckL]] + A1 . E1; \[GothicCapitalC] = 
          InvertRows[A3] . E1; DInv = Inverse[\[GothicCapitalD]]; 
    Omega1 = Inverse[\[GothicCapitalD] - \[GothicCapitalC] . DInv . 
             \[GothicCapitalC]]; Omega2 = DInv . \[GothicCapitalC] . Omega1; 
    W = CovarianceMatrixARMA[theta, phi, n]; 
    S1 = SubMatrix[W, {1, 1}, {\[DoubleStruckL], \[DoubleStruckL]}]; 
    S2 = SubMatrix[W, {\[DoubleStruckL] + 1, 1}, 
           {n - 2*\[DoubleStruckL], \[DoubleStruckL]}]; 
    S3 = SubMatrix[W, {n - \[DoubleStruckL] + 1, 1}, {\[DoubleStruckL], 
            \[DoubleStruckL]}]; S22 = SubMatrix[W, {\[DoubleStruckL] + 1, 
            \[DoubleStruckL] + 1}, {n - 2*\[DoubleStruckL], 
            n - 2*\[DoubleStruckL]}]; T1 = Omega1 . (S1 - A1 . E2) - 
           Omega2 . InvertRows[S3 - A3 . E2]; T2 = S2 . Transpose[Omega1] - 
           InvertRows[S2] . Transpose[Omega2]; 
    T3 = InvertRC[Omega1] . (S3 - A3 . E2) - InvertRows[Omega2] . 
            (S1 - A1 . E2); T22 = S22 - A2 . E1 . Transpose[T2] - 
           InvertRC[A2 . E1 . Transpose[T2]]; 
    icov = BlockMatrix[{{T1, Transpose[T2], Transpose[T3]}, 
            {T2, T22, InvertRC[T2]}, {T3, InvertRC[Transpose[T2]], 
             InvertRC[T1]}}]; 
   If[ifdet, 
          {DetAR[phi]*Det[BlockMatrix[{{\[GothicCapitalD], 
                 \[GothicCapitalC]}, {\[GothicCapitalC], 
                 \[GothicCapitalD]}}]]/DetAR[theta], icov}, 
     icov]]] /; n > 0
 
 
InverseCovarianceMatrixAR[phi_List, n_Integer, (opts___)?OptionQ] := 
    Module[{ifdet, i}, ifdet = Determinant /. {opts} /. 
         Options[InverseCovarianceMatrixAR]; i = SiddiquiMatrix[phi, n]; 
       If[ifdet, {DetAR[phi],i}, i]] /; n > 0

Options[InverseCovarianceMatrixAR] = {Determinant -> False};
 
InverseCovarianceMatrixMA[theta_List, n_Integer, (opts___)?OptionQ] := 
    Module[{q = Length[theta], ifdet, W, A1, A2, A3, A22, theta2, B, 
       \[GothicCapitalD], \[GothicCapitalC], DInv, Omega1, Omega2, M1, M2, 
       M3, M22, icov}, ifdet = Determinant /. {opts} /. 
         Options[InverseCovarianceMatrixMA]; W = CovarianceMatrixARMA[theta, 
         {}, n]; A1 = SubMatrix[W, {1, 1}, {q, q}]; 
       A2 = SubMatrix[W, {q + 1, 1}, {n - 2*q, q}]; 
       A3 = SubMatrix[W, {n - q + 1, 1}, {q, q}]; 
       A22 = SubMatrix[W, {q + 1, q + 1}, {n - 2*q, n - 2*q}]; 
       theta2 = Join[{-1}, theta]; B = Table[0, {q}, {q}]; 
       For[j = 1, j <= q, ++j, For[i = 1, i <= q, ++i, 
         B[[i,j]] = Which[j > i, B[[j,i]], j <= i, 
           Sum[theta2[[1 + q + k - i]]*theta2[[1 + q + k - j]], 
            {k, 1, Min[i, j]}]]]]; B = InvertRC[B]; \[GothicCapitalD] = 
        IdentityMatrix[q] + B . A1; \[GothicCapitalC] = B . InvertRows[A3]; 
       DInv = Inverse[\[GothicCapitalD]]; 
       Omega1 = Inverse[\[GothicCapitalD] - \[GothicCapitalC] . DInv . 
           \[GothicCapitalC]]; Omega2 = DInv . \[GothicCapitalC] . Omega1; 
       M1 = A1 . Omega1 - InvertRows[A3] . Omega2; 
       M2 = A2 . Omega1 - InvertRows[A2] . Omega2; 
       M3 = A3 . Omega1 - InvertRows[A1] . Omega2; 
       M22 = A22 - M2 . B . Transpose[A2] - InvertRows[M2] . B . 
          InvertColumns[Transpose[A2]]; 
       icov = BlockMatrix[{{M1, Transpose[M2], Transpose[M3]}, 
          {M2, M22, InvertRC[M2]}, {M3, InvertRC[Transpose[M2]], 
           InvertRC[M1]}}]; If[ifdet, 
        {Det[BlockMatrix[{{\[GothicCapitalD], \[GothicCapitalC]}, 
             {\[GothicCapitalC], \[GothicCapitalD]}}]]/DetAR[theta], icov}, 
        icov]] /; n > 0
 
Options[InverseCovarianceMatrixMA] = {Determinant -> False};

InvertRC[a_] := InvertColumns[InvertRows[a]];
 
InvertColumns[a_] := Transpose[InvertRows[Transpose[a]]];
 
InvertRows[a_] := With[{n = Length[a]}, Table[a[[n + 1 - j]], {j, n}]];

InformationMatrixARMA[(phi_)?VectorQ, (theta_)?VectorQ, (opts___)?OptionQ] := 
  Module[{vvmatrix, uvmatrix, uumatrix, i, id}, 
   fd = FractionalDifference /. {opts} /. Options[InformationMatrixARMA]; 
    Which[Length[phi] > 0, vvmatrix = 
      Toeplitz[Drop[Drop[tccf[phi, phi], Length[phi] - 1], -1]], True, 
     vvmatrix = {}]; Which[Length[phi] > 0 && Length[theta] > 0, 
     uvmatrix = 
      Table[-tccf[phi, theta][[Length[theta] + ii - jj]], 
       {ii, Length[phi]}, {jj, Length[theta]}], True, vumatrix = {}]; 
    Which[Length[theta] > 0, uumatrix = 
      Toeplitz[Drop[Drop[tccf[theta, theta], Length[theta] - 1], -1]], True, 
     uumatrix = {}]; i = 
     Which[Length[phi] > 0 && Length[theta] > 0, 
      Join[Transpose[Join[vvmatrix, Transpose[uvmatrix]]], 
       Transpose[Join[uvmatrix, uumatrix]]], 
      Length[theta] == 0 && Length[phi] > 0, vvmatrix, 
      Length[phi] == 0 && Length[theta] > 0, uumatrix, True, {{}}]; 
    If[ !fd, Return[i]]; id = 
     Transpose[{Join[FracPhi[phi], FracPhi[theta]]}]; 
    Join[RJoin[i, id], {Flatten[Append[id, {Pi^2/6//N}]]}]];
 
InformationMatrixARMA[(phi_)?VectorQ, (theta_)?VectorQ, (phis_)?VectorQ, 
   (thetas_)?VectorQ, (s_Integer)?NonNegative, (opts___)?OptionQ] := 
  Module[{i1, i2, i12aa, i12ab, i12ba, i12bb, i12, p = Length[phi], 
    q = Length[theta], ps = Length[phis], qs = Length[thetas], i, fd}, 
   fd = FractionalDifference /. {opts} /. Options[InformationMatrixARMA]; 
    i1 = InformationMatrixARMA[phi, theta, FractionalDifference -> False]; 
    i2 = InformationMatrixARMA[phis, thetas, FractionalDifference -> False]; 
    i12aa = If[p > 0 && ps > 0, tccfseasonal[phi, phis, s], {{}}]; 
    i12ab = If[p > 0 && qs > 0, tccfseasonal[phi, thetas, s], {{}}]; 
    i12ba = If[q > 0 && ps > 0, tccfseasonal[theta, phis, s], {{}}]; 
    i12bb = If[q > 0 && qs > 0, tccfseasonal[theta, thetas, s], {{}}]; 
    i12 = Join[RJoin[i12aa, i12ab], RJoin[i12ba, i12bb]]; 
    i = Which[ps == 0 && qs == 0, i1, p == 0 && q == 0, i2, True, 
      Join[RJoin[i1, i12], RJoin[Transpose[i12], i2]]]; If[ !fd, Return[i]]; 
    id = Transpose[{Join[FracPhi[phi], FracPhi[theta], FracPhi[phis, s], 
        FracPhi[thetas, s]]}]; Join[RJoin[i, id], 
     {Flatten[Append[id, {Pi^2/6//N}]]}]];

DetARMA[phi_, theta_] := 
  Which[Length[phi] == 0, DetAR[theta], Length[theta] == 0, DetAR[phi], True, 
   Simplify[DetAR[Drop[-CoefficientList[ARPolynomial[phi, B]*
          ARPolynomial[theta, B], B], 1]]/(DetAR[phi]*DetAR[theta])^2]];
 
DetAR[phi_] := 
  With[{p = Length[phi]}, 
    Inner[Power, (1 - ARToPacf[phi]^2), -Range[p], Times]];

TacvfFDWN[(d_)?NumberQ, MaxLag_Integer /; MaxLag >= 0] := 
  Module[{k = 0}, NestList[(++k - 1 + d)/(k - d)*#1 & , 
     Gamma[1 - 2*d]/Gamma[1 - d]^2, MaxLag]] /; Abs[d] < 0.5;
 
FDtoGLP[d_Real, MaxLag_Integer] := 
  Module[{B}, Drop[CoefficientList[Normal[Series[(1 - B)^(-d), 
       {B, 0, MaxLag}]], B], 1]];

FDtoAR[d_, MaxLag_] := 
  Module[{B}, -Drop[CoefficientList[Normal[Series[(1 - B)^d, 
        {B, 0, MaxLag}]], B], 1]];

 
FracPhi[phi_] := 
  If[Length[phi] == 0, {}, Table[Sum[ARToMACoefficient[phi, k]/(k + i + 1), 
     {k, 0, 250}], {i, 0, Length[phi] - 1}]];
 
FracPhi[phis_, s_] := 
  If[Length[phis] == 0 || s == 0, {}, 
   Table[Sum[ARToMACoefficient[phis, k]/((k + i + 1)*s), {k, 0, 250}], 
    {i, 0, Length[phis] - 1}]];
 
tccfseasonal[phi_, phis_, s_] := 
  Module[{p = Length[phi], ps = Length[phis], c}, 
   If[p == 0 || ps == 0, Return[{}]]; 
    c = tccf[phi, ExpandToSeasonal[phis, s]]; 
    Table[c[[ps*s + ii - jj*s]], {ii, p}, {jj, ps}]];
 
(*duplicate code. remove!
**
*InformationMatrixARMA[(phi_)?VectorQ, (theta_)?VectorQ, (opts___)?OptionQ] := 
*  Module[{vvmatrix, uvmatrix, uumatrix, i, id}, 
*   fd = FractionalDifference /. {opts} /. Options[InformationMatrixARMA]; 
*    Which[Length[phi] > 0, vvmatrix = 
*      Toeplitz[Drop[Drop[tccf[phi, phi], Length[phi] - 1], -1]], True, 
*     vvmatrix = {}]; Which[Length[phi] > 0 && Length[theta] > 0, 
*     uvmatrix = 
*      Table[-tccf[phi, theta][[Length[theta] + ii - jj]], 
*       {ii, Length[phi]}, {jj, Length[theta]}], True, vumatrix = {}]; 
*    Which[Length[theta] > 0, uumatrix = 
*      Toeplitz[Drop[Drop[tccf[theta, theta], Length[theta] - 1], -1]], True, 
*     uumatrix = {}]; i = 
*     Which[Length[phi] > 0 && Length[theta] > 0, 
*      Join[Transpose[Join[vvmatrix, Transpose[uvmatrix]]], 
*       Transpose[Join[uvmatrix, uumatrix]]], 
*      Length[theta] == 0 && Length[phi] > 0, vvmatrix, 
*      Length[phi] == 0 && Length[theta] > 0, uumatrix, True, {{}}]; 
*    If[ !fd, Return[i]]; id = 
*     Transpose[{Join[FracPhi[phi], -FracPhi[theta]]}]; 
*    Join[RJoin[i, id], {Flatten[Append[id, {Pi^2/6//N}]]}]];
* 
*InformationMatrixARMA[(phi_)?VectorQ, (theta_)?VectorQ, (phis_)?VectorQ, 
*   (thetas_)?VectorQ, (s_Integer)?NonNegative, (opts___)?OptionQ] := 
*  Module[{i1, i2, i12aa, i12ab, i12ba, i12bb, i12, p = Length[phi], 
*    q = Length[theta], ps = Length[phis], qs = Length[thetas], i, fd}, 
*   fd = FractionalDifference /. {opts} /. Options[InformationMatrixARMA]; 
*    i1 = InformationMatrixARMA[phi, theta, FractionalDifference -> False]; 
*    i2 = InformationMatrixARMA[phis, thetas, FractionalDifference -> False]; 
*    i12aa = If[p > 0 && ps > 0, tccfseasonal[phi, phis, s], {{}}]; 
*    i12ab = If[p > 0 && qs > 0, tccfseasonal[phi, thetas, s], {{}}]; 
*    i12ba = If[q > 0 && ps > 0, tccfseasonal[theta, phis, s], {{}}]; 
*    i12bb = If[q > 0 && qs > 0, tccfseasonal[theta, thetas, s], {{}}]; 
*    i12 = Join[RJoin[i12aa, i12ab], RJoin[i12ba, i12bb]]; 
*    i = Which[ps == 0 && qs == 0, i1, p == 0 && q == 0, i2, True, 
*      Join[RJoin[i1, i12], RJoin[Transpose[i12], i2]]]; If[ !fd, Return[i]]; 
*    id = Transpose[{Join[FracPhi[phi], -FracPhi[theta], FracPhi[phis, s], 
*        -FracPhi[thetas, s]]}]; 
*    Join[RJoin[i, id], {Flatten[Append[id, {Pi^2/6//N}]]}]];
*)

Toeplitz/: Toeplitz[m_] := 
   Return[Table[m[[1 + Abs[i - j]]], {i, Length[m]}, {j, Length[m]}]];

(* tccf - calculates the theoretical auto & cross correlation functions
 Reference: "Derivation of the theoretical autocorrelation function of
 Autoregressive-Moving Average Time Series, Applied Statistics, 24, 255-6. *) 
tccf/: tccf[phi_, theta_] := 
   Module[{k, rhs, amatrix}, k = Length[phi] + Length[theta]; 
     rhs = Table[0, {k}]; rhs[[1]] = -1; 
     amatrix = Table[getaij[phi, theta, i, j], {i, k}, {j, k}]; 
     Return[LinearSolve[amatrix, rhs]]];

(* getaij - utility function required by tccf*) 
getaij/: getaij[phi_, theta_, (i_Integer)?Positive, (j_Integer)?Positive] := 
   Module[{p, q, imj, ijq}, p = Length[phi]; q = Length[theta]; imj = i - j; 
     ijq = (i + j) - q - 1; Return[Which[i > q && i > j && imj <= q, 
       Rationalize[theta[[imj]]], i > q && imj == 0, -1, 
       i <= q && ijq > 0 && ijq <= p, Rationalize[phi[[ijq]]], 
       i <= q && ijq == 0, -1, True, 0]]];

InformationMatrixFARMA[phi_, theta_, d_, maxlag_] := 
  Module[{I1, J1, J2, J3, I2, I3},
  Which[Length[phi]>0 || Length[theta]>0,
        I1 = InformationMatrixARMA[phi, theta],
        True, I1={}];
  Which[Length[phi]>0,
        J1 = -Table[
             Sum[ARToMACoefficient[phi, k]/(k+i+1), {k,0,maxlag}], 
             {i,0,Length[phi]-1}],
        True,
        {}];
  Which[Length[theta]>0,
        J2 = Table[
             Sum[ARToMACoefficient[theta, k]/(k+i+1), {k,0,maxlag}], 
             {i,0,Length[theta]-1}],
        True,
        {}];
  J3 = Pi^2/6;
  I2 = Join[I1, {Join[J1,J2]}];
  I3 = Transpose[Join[Transpose[I2], {Join[J1,J2,{J3}]}]];
  Return[I3]
  ];



SiddiquiMatrix[phi_List, n_Integer/;n>0] := 
  Module[{largephi, p, A}, 
     largephi = Join[{-1}, phi, Array[0&,{n-Length[phi]}]]; 
     p    = Length[largephi]-1; 
     A    = Table[0,{p},{p}];
     For[j = 1, j <= n, ++j, 
     For[i = 1, i <= n, ++i, 
           A[[i, j]] = 
                     Which[j >  i,  A[[j, i]], 
                           j <= i, 
            Sum[largephi[[1+i-k]]*largephi[[1+j-k]]-
                largephi[[1+p+k-i]]*largephi[[1+p+k-j]]
	     ,{k,1,Min[i,j]}]
	   ]
     ]];
	 A
  ];

SiddiquiMatrix[phi_Symbol, n_Integer/;n>0] := SiddiquiMatrix[{phi}, n];

SiddiquiMatrix[phi_Number, n_Integer/;n>0] := SiddiquiMatrix[{phi}, n];

(*=========================ACF==============================*)
PlotRacf[res_, phi_, theta_, (opts___)?OptionQ] := 
  Module[{maxlag, acfrange, frameticks, r, blim, sbound, squared}, 
   maxlag = MaxLag /. {opts} /. Options[PlotRacf]; 
    sbound = SimultaneousBound /. {opts} /. Options[PlotRacf]; 
    blim = If[sbound, Quantile[NormalDistribution[0, 1], 0.975^(1/maxlag)], 
      1.96]; squared = 
     SquaredResiduals /. {opts} /. Options[PlotRacf]; 
    r = If[squared, Sacf[res*res, MaxLag -> maxlag, MeanZero -> False], 
      Sacf[res, MaxLag -> maxlag, MeanZero -> True]]; 
    sd = If[squared, Table[1/Sqrt[Length[res]], {maxlag}], 
      Sqrt[Transpose[VarianceRACFARMA[phi, theta, maxlag, Length[res]], 
        {1, 1}]]]; acfrange = 
     {-1, 1}*Ceiling[10*1.06*Max[Abs[r], blim/Sqrt[Length[res]]]]/10.; 
    frameticks = 
     {Range[1, Length[r], Max[1, Round[Length[r]/10.]]], 
      Range[acfrange[[1]], acfrange[[2]], 0.1]}//Chop;
    Show[Graphics[{{Line[{{1, 0}, {Length[r], 0}}]}, 
       {RGBColor[0, 0, 1], Thickness[0.01], 
        Line /@ ({{#1, 0}, {#1, r[[#1]]}} & ) /@ Range[Length[r]]}, 
       {RGBColor[1, 0, 0], Thickness[0.01], 
        Line[Transpose[{Range[Length[r]], blim*sd}]], 
        Line[Transpose[{Range[Length[r]], -blim*sd}]]}}, Frame -> True, 
      Axes -> False, AxesOrigin -> {1, 0}, FrameLabel -> {"Lag", "Acf"}, 
      PlotRange -> 
       With[{eps = 0.05*maxlag}, {{1 - eps, maxlag + eps}, acfrange}], 
      FrameTicks -> frameticks]]];

VarianceRacfARMA[phi_, theta_, MaxLag_, n_] := 
  Module[{X = {}, x, y, p = Length[phi], q = Length[theta]}, 
   If[p > 0, x = -Prepend[ARInvert[phi, MaxLag - 1], 1]]; 
    If[q > 0, y = Prepend[ARInvert[theta, MaxLag - 1], 1]]; 
    Do[AppendTo[X, Join[Table[0, {i}], Drop[x, -i]]], {i, 0, p - 1}]; 
    Do[AppendTo[X, Join[Table[0, {i}], Drop[y, -i]]], {i, 0, q - 1}]; 
    Simplify[(IdentityMatrix[MaxLag] - 
    Transpose[X] . Inverse[InformationMatrixARMA[phi, theta]] . X)/n]];
 
Sacf[(z_)?VectorQ, (opts___)?OptionQ] := 
  Module[{ifmean0, maxlag, x, X, c}, 
   ifmean0 = MeanZero /. {opts} /. Options[Sacf]; 
    maxlag = MaxLag /. {opts} /. Options[Sacf]; 
    If[maxlag <= 0, Print["MaxLag<0"]; Return[{}]]; 
    X = {x0 = x = z - If[ifmean0, 0, Mean[z]]}; 
    If[Plus @@ Abs[x] == 0, Print["constant series!"]; Return[{}]]; 
    Do[AppendTo[X, x = Drop[Prepend[x, 0], -1]], {maxlag}]; c = X . x0; 
    Drop[c, 1]/First[c]];
 
Sacvf[(z_)?VectorQ, (opts___)?OptionQ] := 
  Module[{ifmean0, maxlag, x, X}, 
   ifmean0 = MeanZero /. {opts} /. Options[Sacvf]; 
    maxlag = MaxLag /. {opts} /. Options[Sacvf]; 
    If[maxlag <= 0, Print["MaxLag<0"]; Return[{}]]; 
    X = {x0 = x = z - If[ifmean0, 0, Mean[z]]}; 
    If[Plus @@ Abs[x] == 0, Print["constant series!"]; Return[{}]]; 
    Do[AppendTo[X, x = Drop[Prepend[x, 0], -1]], {maxlag}]; X . x0/Length[z]
    ];

SacfPlot[z_,opts___?OptionQ] :=
   Module[{ifmean0,maxlag,cilim,frameticks,r},
		ifmean0=MeanZero/.{opts}/.Options[SacfPlot];
		maxlag=MaxLag/.{opts}/.Options[SacfPlot];
		blag=BartlettLargeLag/.{opts}/.Options[SacfPlot];
    sbound=Bound/.{opts}/.Options[SacfPlot];
    plab=PlotLabel/.{opts}/.Options[SacfPlot];
			r=Prepend[Sacf[z,MaxLag->maxlag,MeanZero->ifmean0],1];
			cilim=1.96*If[blag>0,Apply[Plus,Take[r,blag]^2],1]/Sqrt[Length[z]];
		slim=Quantile[NormalDistribution[0,1],0.975^(1/maxlag)]*cilim/1.96;
		frameticks={Range[0,Length[r],Max[1,Round[Length[r]/10.0]]], 
        Range[-1,1,.5]};
   Show[Graphics[
    {{RGBColor[0,0,1], Thickness[0.01],
     Map[Line, 
         Map[{{#1 - 1, 0}, {#1 - 1, 
                      r\[LeftDoubleBracket]#1\[RightDoubleBracket]}} & , 
                Range[Length[r]]]
        ]
    },
         {Line[{{0,0},{maxlag,0}}]},
    Which[
MatchQ[sbound,"Individual"],{RGBColor[1,0,0],Thickness[0.01],
              Line[{{blag+1,cilim},{maxlag,cilim}}],
              Line[{{blag+1,-cilim},{maxlag,-cilim}}]},
MatchQ[sbound,"Simultaneous"],{RGBColor[0,1,0],
              Thickness[0.01],Line[{{blag+1,slim},{maxlag,slim}}],
              Line[{{blag+1,-slim},{maxlag,-slim}}]},
True,{}]},
                Frame->True,
	PlotLabel->plab,
				Axes->False,
                AxesOrigin->{0,0}, 
                FrameLabel->{"Lag", "Acf"},
                PlotRange->{-1,1},
                FrameTicks->frameticks
            ]]
	];

(*=========================PACF==============================*)
PacfDL[(c_)?NumberVectorQ, (opts___)?OptionQ] := 
  Module[{a, phik, vk, phikm1, vkm1, L = Length[c] - 1, 
    d = Drop[c, 1], pi, ifar}, 
     ifar = LinearPredictor /. {opts} /. Options[PacfDL]; 
     If[L==0,phik={ };vk=c[[1]];pi={ },
     phik = {c[[2]]/c[[1]]};
     pi = Table[0, {L}]; 
    pi[[1]] = First[phik]; vk = c[[1]]*(1 - pi[[1]]^2);]; 
    If[L>1,Do[vkm1 = vk; phikm1 = phik; a = DLa[phikm1, d, k]/vk; 
      phik = DLphik[phikm1, a]; vk = vkm1*(1 - a^2); pi[[k]] = a, 
     {k, 2, L}]]; 
    If[ifar, {ARCoefficients -> phik, ResidualVariance -> vk, Pacf -> pi}, 
             pi]];
 
PacfDL[c_List, (opts___)?OptionQ] := 
  Module[{a, phik = { }, vk, phikm1, vkm1, L = Length[c] - 1, 
    d = Drop[c, 1], pi, ifar}, 
    ifar = LinearPredictor /. {opts} /. Options[PacfDL]; 
    If[L==0,phik={};vk=c[[1]];pi={},
    phik = {c[[2]]/c[[1]]};
    pi = Table[0, {L}]; 
    pi[[1]] = First[phik]; vk = c[[1]]*(1 - pi[[1]]^2);];
    If[L>1,Do[vkm1 = vk; phikm1 = phik; 
      a = Prepend[-phikm1, 1] . Reverse[Take[d, k]]/vk; 
      phik = Append[phikm1 - a*Reverse[phikm1], a]; vk = vkm1*(1 - a^2); 
      pi[[k]] = a, {k, 2, L}]]; If[ifar, 
     {ARCoefficients -> phik, ResidualVariance -> vk, Pacf -> pi}, pi]];
 
PacfBayes[z_?NumberVectorQ, MaxLag_Integer, (opts___)?OptionQ] := 
  Module[{}, ifar = LinearPredictor /. {opts} /. Options[PacfBayes]; 
    pacf = PacfBurg[z, MaxLag]; 
    aic = ICFromPacf[Length[z], p, Table -> True, Criterion -> "AIC"]; 
    Lp = Exp[-0.5*(aic - Min[aic])]; 
    prior = 1./With[{k = Range[0, MaxLag]}, (k + 1)*(k + 2)]; 
    posterior = Lp*prior/Plus @@ (Lp*prior); 
    wts = 1 - Drop[Drop[FoldList[Plus, 0, posterior], 1], -1]; 
    pacfbayes = pacf*wts; If[ifar, 
     {ARCoefficients -> PacfToAR[pacfbayes], 
      ResidualVariance -> VarianceMLE[z]*Times @@ (1 - pacfbayes^2), 
      Pacf -> pacfbayes}, pacfbayes]] /; MaxLag>0 ;

PacfPlot[(z_)?NumberVectorQ, (opts___)?OptionQ] := 
  Module[{slim,sbound,method, ifmean0, maxlag, cilim, frameticks, r}, 
   method = Method /. {opts} /. Options[PacfPlot]; 
    ifmean0 = MeanZero /. {opts} /. Options[PacfPlot]; 
    maxlag = MaxLag /. {opts} /. Options[PacfPlot]; 
    sbound = SimultaneousBound /. {opts} /. Options[SacfPlot]; 
    r = Prepend[Which[MatchQ[method, "Burg"], 
       PacfBurg[z, maxlag, MeanZero -> ifmean0], True, 
       PacfDL[Sacvf[z, MaxLag -> maxlag, MeanZero -> ifmean0]]], 1]; 
    cilim = 1.96/Sqrt[Length[z]]; 
    slim = Quantile[NormalDistribution[0, 1], 0.975^(1/maxlag)]*
      cilim/1.96; frameticks = 
     {Range[0, Length[r], Max[1, Round[Length[r]/10.]]], Range[-1, 1, 0.5]}; 
    Show[Graphics[{{RGBColor[0, 0, 1], Thickness[0.01], 
        Line /@ ({{#1 - 1, 0}, {#1 - 1, r[[#1]]}} & ) /@ Range[Length[r]]}, 
       If[ !sbound, {RGBColor[1, 0, 0], Thickness[0.01], 
         Line[{{1, cilim}, {maxlag, cilim}}], 
         Line[{{1, -cilim}, {maxlag, -cilim}}]}, 
        {RGBColor[1, 1, 0], Thickness[0.01], 
         Line[{{1, slim}, {maxlag, slim}}], 
         Line[{{1, -slim}, {maxlag, -slim}}]}, {}]}, Frame -> True, 
      Axes -> False, AxesOrigin -> {0, 0}, FrameLabel -> {"Lag", "Acf"}, 
      PlotRange -> {-1, 1}, FrameTicks -> frameticks, 
      GridLines -> {{}, Range[-1, 1, 0.5]}]]];

ARToPacf[phi_?VectorQ] := 
  Module[{pi, a, phikp1, phik = phi, L = Length[phi]}, 
   pi = Table[0, {L}]; Do[pi[[L + 1 - k]] = a = Last[phik]; 
      phikp1 = Drop[phik, -1]; phik = (phikp1 + a*Reverse[phikp1])/(1 - a^2)\
, {k, 1, L}]; pi];
 
PacfToAR[pi_?VectorQ] := 
  Module[{L = Length[pi], phik, phikm1}, 
   phik = {pi[[1]]}; Do[phikm1 = phik; 
      phik = Append[phikm1 - pi[[k]]*Reverse[phikm1], pi[[k]]], {k, 2, L}]; 
    phik];
 
ICFromPacf[n_Integer, (pacf_)?NumberVectorQ, (opts___)?OptionQ] := 
  Module[{tableQ,whichic, neg2LL, plotQ}, whichic = 
      Criterion /. {opts} /. Options[ICFromPacf]; 
     tableQ = Plot /. {opts} /. Options[ICFromPacf]; 
     plotQ = Plot /. {opts} /. Options[ICFromPacf]; 
     neg2LL = n*Log[FoldList[#1*(1 - #2^2) & , 1, pacf]]; 
     ic = Which[StringMatchQ[whichic, "AIC"], 
       neg2LL + 2*Range[0, Length[pacf]], True, 
       neg2LL + Log[n]*Range[0, Length[pacf]]]; 
     Which[tableQ,ic,plotQ, ListPlot[Transpose[{Range[0, Length[pacf]], ic}], 
       Axes -> False, Frame -> True, PlotRange -> All, 
       PlotStyle -> {AbsolutePointSize[6]}], 
       True, IndexOfMinimum[ic]-1]]/; n > 0;
 
PacfBurg[(z_)?NumberVectorQ, MaxLag_Integer, (opts___)?OptionQ] := 
  Module[{meanzero, ific, lp, zc, fk1, bk1, fk0, bk0, 
          pi = Table[0, {MaxLag}], L}, 
    meanzero = MeanZero /. {opts} /. Options[PacfBurg]; 
     ific = ModelSelection /. {opts} /. Options[PacfBurg]; 
     lp = LinearPredictor /. {opts} /. Options[PacfBurg]; 
     fk1 = bk1 = zc = If[meanzero, z, z - Mean[z]]; 
     Do[pi[[k]] = 
        With[{f = Drop[fk1, 1], b = Drop[bk1, -1]}, 
         2*f . b/Plus @@ (Join[f, b]^2)]; fk0 = fk1; bk0 = bk1; 
       fk1 = Drop[fk0, 1] - pi[[k]]*Drop[bk0, -1]; 
       bk1 = Drop[bk0, -1] - pi[[k]]*Drop[fk0, 1], {k, MaxLag}]; 
     L = If[ific, ICFromPacf[pi, Null], MaxLag]; 
     If[ !lp, pi, If[L == 0, phik = {}; vk = Mean[zc^2], 
        phik = {pi[[1]]}; vk = Mean[zc^2]*(1 - pi[[1]]^2)]; 
       If[L > 1, Do[vkm1 = vk; phikm1 = phik; 
          phik = Append[phikm1 - pi[[k]]*Reverse[phikm1], pi[[k]]]; 
          vk = vkm1*(1 - pi[[k]]^2), {k, 2, L}]]; 
       {ARCoefficients -> phik, ResidualVariance -> vk, Pacf -> pi}]] /; 
   MaxLag > 0;

(*=========================MA(1) Estimation==============================*)

MA1Fit[(zz_)?NumberVectorQ, (opts___)?OptionQ] := 
    Module[{ETOL = 0.0001, NULLVEC = {{NULL, NULL}, {NULL, NULL}, 
        {NULL, NULL}}, ls, \[Theta], maxls, i, a, fa, b, fb, c, fc, d, fd, 
      test1, test2, fx, x, z = N[zz]}, 
     z = If[MeanZero /. {opts} /. Options[MA1Fit], z, z - Mean[z]]; 
      ls = Table[{\[Theta], logLMA1F[\[Theta], z]}, {\[Theta], -1, 1, 0.1}]; 
      {maxls} = Select[ls, #1[[2]] == Max[Transpose[ls][[2]]] & ]; 
      i = Floor[maxls[[1]]*10 + 11]; {{a, fa}, {b, fb}, {c, fc}} = 
       Which[i == 1, d = -0.999; fd = logLMA1F[d, z]; If[ls[[1]][[2]] > fd, 
          NULLVEC, {ls[[1]], {d, fd}, ls[[2]]}], i == 21, 
        d = 0.999; fd = logLMA1F[d, z]; If[ls[[21]][[2]] > fd, NULLVEC, 
          {ls[[20]], {d, fd}, ls[[21]]}], i > 1 && i < 21, 
        {{a, fa}, {b, fb}, {c, fc}} = ls[[{i - 1, i, i + 1}]]]; 
      (If[a === NULL, Which[i == 1, b = -1., i == 21, b = 1.], 
        While[(test1 = c - a > ETOL) || (test2 = fb - Max[fa, fc]) > ETOL, 
         fx = logLMA1F[x = b + 0.5*(((b - a)^2*(fb - fc) - (b - c)^2*
                  (fb - fa))/((b - a)*(fb - fc) - (b - c)*(fb - fa))), z]; 
          {{a, fa}, {b, fb}, {c, fc}} = If[fx > fb, Which[x < b, 
             {{a, fa}, {x, fx}, {b, fb}}, x > b, {{b, fb}, {x, fx}, 
              {c, fc}}], If[b - a > 0.5*(c - a), d = b + 0.38196*(a - b); 
              fd = logLMA1F[d, z]; If[fd > fb, {{a, fa}, {d, fd}, {b, fb}}, {
                {d, fd}, {b, fb}, {c, fc}}], d = b + 0.38196*(c - b); 
              fd = logLMA1F[d, z]; If[fd > fb, {{b, fb}, {d, fd}, {c, fc}}, {
                {a, fa}, {b, fb}, {d, fd}}]]]]]; b)];

logLMA1F[\[Theta]_, z_] := Module[{n = Length[z]}, 
     detma = DetMA[\[Theta], n]; Lz = GetLz[\[Theta], z]; 
      u = Getu0[\[Theta], Lz, detma]; S = GetSumSq[\[Theta], Lz, u, n]; 
      (-(1/2))*Log[detma] - (n/2)*Log[S/n]];
 
GetLz=Compile[{{\[Theta],_Real},{z,_Real,1}},FoldList[(#1 \[Theta]+#2)&,0.0,z]];

Getu0=Compile[{{\[Theta],_Real},{Lz,_Real,
          1},{detma,_Real}},-Fold[#1 \[Theta]+#2&,0.0,Reverse[Lz]]/detma];

DetMA=Compile[{{\[Theta],_Real},{n,_Integer}},
      Apply[Plus,NestList[#1 \[Theta]&,1,n]^2]];

GetSumSq=Compile[{{\[Theta],_Real},{Lz,_Real,1},{u,_Real},{n,_Integer}},
      Apply[Plus,(Lz+NestList[#1 \[Theta]&,1,n] u)^2]];

 
(*=========================AR(1) Estimation==============================*)

getD11 = Compile[{{z, _Real, 1}}, z.z];

getD22 = Compile[{{z, _Real, 1}, {d11, _Real}},
      d11 - 
        With[{y = z\[LeftDoubleBracket]{1, Length[z]}\[RightDoubleBracket]}, 
          y.y]];

getD12 = Compile[{{z, _Real, 1}}, Drop[z, {1}].Drop[z, {-1}]];

getmean = 
    Compile[{{x, _Real, 
          1}, {\[Rho], _Real}}, (Plus @@ 
              x - \[Rho] Plus @@ Drop[Drop[x, {1}], -1])/(Length[x] + (2 - Length[x]) \[Rho])];


logLikAR1 := Compile[{{\[Phi], _Real},
      {D11, _Real}, {D12, _Real}, 
      {D22, _Real}, {n, _Integer}}, 
      0.5*(Log[1 - \[Phi]^2] - n*Log[2.*Pi] - n*Log[D11 - 2*\[Phi]*D12 + \[Phi]^2*D22])];
 
logLikRW := Compile[{{z, _Real, 1}, {n, _Integer}}, 
     With[{w = Drop[z, 1] - Drop[z, -1], n = Length[z]}, 
      -0.5*n*(Log[2.0*Pi] + Log[w . w])]];

get\[Phi] = Compile[{{D11, _Complex}, {D12, _Complex}, {D22, _Complex}},
Chop[D12/(3*
                D22) + (((1 - I*Sqrt[3])*(-D12^2 - 3*D22 - 
                        3*D22^2))/(3*2^(2/3)*
                    D22*(2*D12^3 + 9*D12*D22 - 18*D12*D22^2 + 
                          Sqrt[4*(-D12^2 - 3*D22 - 3*D22^2)^3 + (2*D12^3 - 
                                    9*D12*(-1 - D22)*D22 - 
                                    27*D12*D22^2)^2])^(1/3)) - 
              1/(6*2^(1/3)*D22)*(1 + 
                    I*Sqrt[3])*(2*D12^3 + 9*D12*D22 - 18*D12*D22^2 + 
                      Sqrt[4*(-D12^2 - 3*D22 - 3*D22^2)^3 + (2*D12^3 - 
                                9*D12*(-1 - D22)*D22 - 27*D12*D22^2)^2])^(1/
                      3))]
      ];

MarkovModel[(z_)?NumberVectorQ, (opts___)?OptionQ] := 
    Module[{d11, d12, d22, \[Rho] = 0., \[Rho]0, etol, abserr = 100., t, rl, 
      aloglar1, aloglrw, phiopt, indzero, ETOL, indav}, 
     indzero = MeanZero /. {opts} /. Options[MarkovModel]; 
      ETOL = ErrorTolerance /. {opts} /. Options[MarkovModel]; 
      indav = MeanAverage /. {opts} /. Options[MarkovModel]; 
      t = If[indzero, 0, Mean[z]]; d11 = getD11[z - t]; d12 = getD12[z - t]; 
      d22 = getD22[z - t, d11]; \[Rho] = Chop[get\[Phi][d11, d12, d22]]; 
      If[ !indav &&  !indzero, While[abserr > ETOL, 
        t = If[indzero, 0, getmean[z, \[Rho]0 = \[Rho]]]; 
         d11 = getD11[z - t]; d12 = getD12[z - t]; d22 = getD22[z - t, d11]; 
         \[Rho] = Chop[get\[Phi][d11, d12, d22]]; 
         abserr = Abs[\[Rho] - \[Rho]0]/(Abs[\[Rho]] + ETOL); ]]; 
      aloglar1 = logLikAR1[\[Rho], d11, d12, d22, Length[z]]/Length[z]; 
      aloglrw = logLikRW[z - t, Length[z]]/(Length[z] - 1); 
      phiopt = If[aloglar1 > aloglrw, \[Rho], 1]; 
      rl = E^(aloglrw - aloglar1); {\[Phi]hat -> \[Rho], \[Mu]hat -> t, 
       avloglikelihoodar1 -> aloglar1, avloglikelihoodrw -> aloglrw, 
       \[Phi]opt -> phiopt, RLunitroot -> rl}];

MarkovFit[z_List, (opts___)?OptionQ] := 
    Module[{w, phiar1, phiopt, fullreport, llrw, llar1, tmp}, 
     w = If[MeanZero /. {opts} /. Options[MarkovFit], z, z - Mean[z]]; 
      fullreport = FullReport /. {opts} /. Options[MarkovFit]; 
      phiar1 = AR1Fit[w, MeanZero -> True]; 
      phiopt = If[(llar1 = AverageLoglikelihoodMarkov[w, phiar1]) < 
         (llrw = AverageLoglikelihoodRandomWalk[z]), 1, phiar1]; 
      If[ !fullreport, phiopt, {phi -> phiopt, loglikelihoodar1 -> llar1, 
        loglikelihoodrw -> llrw, phimle -> phiar1}]];

AverageLoglikelihoodMarkov[z_List, (phi_)?NumberQ] := 
    (LoglikelihoodGP[TacvfARMA[{phi}, {}, Length[z] - 1], z] - 
      0.5*Length[z]*(1 + Log[2*Pi]))/Length[z];
 
AverageLoglikelihoodRandomWalk[z_] := 
    (LoglikelihoodGP[Join[{1}, Table[0, {Length[z] - 1}]], TSDifference[z]] - 
      0.5*(Length[z] - 1)*(1 + Log[2*Pi]))/(Length[z] - 1);
 
AR1Fit[u_List, (opts___)?OptionQ] := Module[{a, b, c, n, t, r, s}, 
      z = If[MeanZero /. {opts} /. Options[AR1Fit], u, u - Mean[u]]; 
      n = Length[z]; a = z . z; 
      b = Drop[z, {1}] . Drop[z, {-1}]; c = a - z[[1]]^2 - z[[n]]^2; 
      Solve[(n - 1)/n*c*t^3 - (n - 2)/n*b*t^2 - (c + a/n)*t + b == 0, t][[
          2]][[1]][[2]]];

OrdinaryLeastSquaresARFit[(z_)?NumberVectorQ, p_Integer /; p > 0] := 
  Module[{}, xs = Table[Subscript[x, j], {j, 1, p}]; 
    ans = Fit[Reverse /@ 
       Drop[Transpose[Reverse[NestList[RotateLeft[#1] & , z, p]]], -p], xs, 
      xs]; Reverse[Table[Coefficient[ans, xs[[j]]], {j, p}]]];
 
ResidualsAR[(z_)?NumberVectorQ, (phihat_)?NumberVectorQ] := 
  With[{p = Length[phihat]}, Drop[Transpose[Reverse[NestList[
        RotateLeft[#1] & , z, p]]], -p] . Prepend[-phihat, 1]];

(*=========================AR Estimation==============================*)
YuleWalkerEstimates[z_?NumberVectorQ, p_Integer/;p>0, opts___?OptionQ]:=
Module[{ifmean0},
ifmean0=MeanZero/.{opts}/.Options[YuleWalkerEstimates];	
PacfDL[Sacvf[z,MaxLag->p,MeanZero->ifmean0],LinearPredictor->True]];

LoglikelihoodAR[phi_, z_, opts___?OptionQ] := 
    Module[{ifmean0, indlags, n = Length[z], p = Length[phi]}, 
     ifmean0 = MeanZero /. {opts} /. Options[LoglikelihoodAR]; 
      w = z - If[ifmean0, 0, Mean[z]]; 
      S = With[{phi2 = Prepend[phi, -1]}, 
           phi2 . ChampernowneD[w, p, MeanZero -> ifmean0] . phi2]; 
      0.5*Log[Det[SiddiquiMatrix[phi, p]]] - 0.5*n*Log[S/n]];
 
ChampernowneD[z_, p_, opts___?OptionQ] := 
    Module[{n = Length[z], c, ifmean0, X, x0, x, \[GothicCapitalE], A}, 
     If[Plus @@ Abs[x] == 0, Print["constant series!"]; Return[{}]]; 
      ifmean0 = MeanZero /. {opts} /. Options[ChampernowneD]; 
      X = {x0 = x = z - If[ifmean0, 0, Mean[z]]}; 
      Do[AppendTo[X, x = Drop[Prepend[x, 0], -1]], {p}]; c = X . x0; 
      \[GothicCapitalE] = Array[#1*0 & , {p + 1, p + 1}]; 
      A = Table[c[[Abs[i - j] + 1]], {i, p + 1}, {j, p + 1}]; 
      Do[\[GothicCapitalE][[i + 1,j + 1]] = \[GothicCapitalE][[i,j]] + 
         z[[i]]*z[[j]] + z[[n + 1 - i]]*z[[n + 1 - j]], {j, p}, {i, j}]; 
      Do[\[GothicCapitalE][[i + 1,j + 1]] = \[GothicCapitalE][[j + 1,i + 1]], 
       {i, p}, {j, i - 1}]; A - \[GothicCapitalE]];

ChampernowneLS[z_, p_, opts___?OptionQ] := 
    Module[{\[GothicCapitalD], ifmean0, indlags, constrainedQ, p2, i = 0}, 
     If[Plus @@ Abs[x] == 0, Print["constant series!"]; Return[{}]]; 
      ifmean0 = MeanZero /. {opts} /. Options[ChampernowneLS]; 
      indlags = Lags /. {opts} /. Options[ChampernowneLS]; 
      constrainedQ = Constrained /. {opts} /. Options[ChampernowneLS]; 
      p2 = If[VectorQ[indlags], Max @@ indlags, p]; \[GothicCapitalD] = 
       If[ifmean0, ChampernowneD[z, p2, MeanZero -> True], 
        ChampernowneD[z, p2, MeanZero -> False]]; \[GothicCapitalD] = 
       If[VectorQ[indlags], lags = Prepend[indlags + 1, 1]; 
         Transpose[Transpose[\[GothicCapitalD][[lags]]][[lags]]], 
        \[GothicCapitalD]]; beta = LinearSolve[Transpose[
         Drop[Transpose[Drop[\[GothicCapitalD], 1]], 1]], 
        Drop[First[\[GothicCapitalD]], 1]]; 
      If[constrainedQ &&  !InvertibleQ[beta], 
       Message[ChampernowneLS::inadmissible]; While[ !InvertibleQ[beta], 
         beta = beta*0.99]; ]; beta];

(*=========================ARMA Estimatoin==============================*)
LoglikelihoodARMA[phi_List, theta_List, z_List, (opts___)?OptionQ] := 
    Module[{g, ifmean0, y, n = Length[z], fd}, 
     ifmean0 = MeanZero /. {opts} /. Options[LoglikelihoodARMA]; 
     fd = FractionalDifference /. {opts} /. Options[LoglikelihoodARMA]; 
     y = z - If[ifmean0, 0, Mean[z]]; 
     LoglikelihoodGP[TacvfARMA[phi, theta, Length[z]-1, FractionalDifference->fd], z]
     ];
 
ResidualsGP[c_?NumberVectorQ, z_?NumberVectorQ, opts___?OptionQ] := 
    Module[{ifmean0, a, phik, phi1, vk, phikm1, vkm1, n, d, S, res}, 
     ifmean0 = MeanZero /. {opts} /. Options[LoglikelihoodGP];
     w = z - If[ifmean0, 0, Mean[z]];
     n = Length[w]; 
     res = Array[0&, n];
     vs = phis = {}; phik = {phi1 = c[[2]]/c[[1]]}; 
     d = Drop[c, 1]; 
     d = With[{m = Length[c]}, 
         Which[n == m, d, 
               n > m,  Join[d, Table[0., {n - m}]], 
               n < m,  Take[d, n - 1]]]; 
     vk = c[[1]]*(1 - phi1^2); 
     res[[1]] = w[[1]]/Sqrt[c[[1]]];  
     res[[2]] = (w[[2]] - phi1*w[[1]])/Sqrt[vk]; 
     g = c[[1]]*vk; 
     Do[vkm1 = vk; phikm1 = phik; phis = Join[phis, {phik}]; 
        a = DLa[phikm1, d, k - 1]/vk; phik = DLphik[phikm1, a]; 
        vk = vkm1*(1 - a^2); g *= vk; vs = Join[vs, vk]; 
        (res[[k]] = DLe[w, phik]/Sqrt[vk]; ), 
       {k, 3, n}]; 
     res];
 
LoglikelihoodGP[c_?NumberVectorQ, z_?NumberVectorQ, opts___?OptionQ] := 
    Module[{ifmean0, a, phik, phi1, vk, phikm1, vkm1, n, d, S}, 
     ifmean0 = MeanZero /. {opts} /. Options[LoglikelihoodGP];
     w = z - If[ifmean0, 0, Mean[z]];
     n = Length[w]; vs = phis = {}; phik = {phi1 = c[[2]]/c[[1]]}; 
     d = Drop[c, 1]; 
     d = With[{m = Length[c]}, 
         Which[n == m, d, 
               n > m,  Join[d, Table[0., {n - m}]], 
               n < m,  Take[d, n - 1]]]; 
     vk = c[[1]]*(1 - phi1^2); 
     S = w[[1]]^2/c[[1]] + 
        (w[[2]] - phi1*w[[1]])^2/vk; 
     g = c[[1]]*vk; 
     Do[vkm1 = vk; phikm1 = phik; phis = Join[phis, {phik}]; 
        a = DLa[phikm1, d, k - 1]/vk; phik = DLphik[phikm1, a]; 
        vk = vkm1*(1 - a^2); g *= vk; vs = Join[vs, vk]; 
        (S += DLe[w, phik]^2/vk; ), 
       {k, 3, n}]; 
     -2^(-1)*Log[g] - n/2*Log[S/n]];

LoglikelihoodGP[c_List, z_List, opts___?OptionQ] := 
      Module[{ifmean0, a, phik = {}, vk, phikm1, vkm1,d = Drop[c, 1], n = Length[z]}, 
      ifmean0 = MeanZero /. {opts} /. Options[LoglikelihoodGP];
      w = z - If[ifmean0, 0, Mean[z]];
      vk = c[[1]]; phik = {c[[2]]/c[[1]]}; 
      vk = c[[1]]*(1 - First[phik]^2); 
      S = w[[1]]^2/c[[1]] + (w[[2]] - First[phik]*w[[1]])^2/vk; 
      g = c[[1]]*vk; 
      Do[vkm1 = vk; phikm1 = phik; 
         a = Prepend[-phikm1, 1] . Reverse[Take[d, k]]/vk; 
         phik = Append[phikm1 - a*Reverse[phikm1], a]; 
         vk = vkm1*(1 - a^2); g *= vk; 
         S += (Prepend[-phik, 1] . Reverse[Take[w, k + 1]])^2/vk, 
       {k, 2, n - 1}]; 
      -2^(-1)*Log[g] - n/2*Log[S/n]];
 

DLa = Compile[{{phikm1,_Real,1},{d,_Real,1},{k, _Integer}},
    Prepend[-phikm1,1.0].Reverse[Take[d,k]]];

DLphik = Compile[{{phikm1,_Real,1},{a,_Real}},
    Append[phikm1-a*Reverse[phikm1],a]];

DLe = Compile[{{z,_Real,1},{phiik,_Real,1}},
    Prepend[-phiik,1.0].Reverse[Take[z,1+Length[phiik]]]];

DLz = Compile[{{z,_Real,1},{phiik,_Real,1},{vk, _Real}},
    Prepend[phiik,Sqrt[vk]].Reverse[Take[z,1+Length[phiik]]]];

(*=========================Panel Time Series==============================*)
FitPanelTimeSeriesAR1[z_?PanelTimeSeriesQ] := 
 Module[{ans, LogLNull, CS, phi1, z2}, 
     z2=Select[z, Not[Apply[Equal, #]] &];
     If[Length[z2]<Length[z], Message[FitPanelTimeSeriesAR1::constant,Length[z2]]];
     ans = FindMinimum[
            -LoglikelihoodPanelTimeSeries[
               TacvfARMA[{phi1}, {}, Max @@ Length /@ z2 - 1], 
               z2], 
            {phi1, {-0.1, 0.1}, -0.999, 0.999}
           ]; 
      LogLNull = LoglikelihoodPanelTimeSeries[
                   TacvfARMA[{0.}, {}, Max @@ Length /@ z2 - 1], z2]; 
      CS = -2*(ans[[1]] + LogLNull); 
      Flatten[{
               phi -> phi1/.ans[[2]], 
               MaximizedLogLikelihood -> -ans[[1]], 
               ChiSqNull -> CS, 
               SignificanceLevel -> 1. - CDF[ChiSquareDistribution[1], CS]}]
 ];
 
LoglikelihoodPanelTimeSeries[c_?NumberVectorQ, z_?PanelTimeSeriesQ] := 
    Plus @@ Table[
      LoglikelihoodGP[c, z[[i]], MeanZero -> True], 
      {i, Length[z]}];

(*=========================Intervention Analysis==========================*)
 
IADesign[(n_)?NumberQ, (T_)?NumberQ, (g_)?NumberQ, (opts___)?OptionQ] := 
    Module[{h, x, type, mean0, dif, dlta}, 
     type = IAType /. {opts} /. Options[IADesign]; 
      dif = Differencing /. {opts} /. Options[IADesign]; 
      dlta = Delta1 /. {opts} /. Options[IADesign]; IADesign[n, {T}, {g}, 
       IAType -> type, Differencing -> dif, Delta1 -> dlta]];
 
IADesign[(n_)?NumberQ, (T_)?NumberVectorQ, (g_)?NumberVectorQ, 
     (opts___)?OptionQ] := Module[{h, x, y, type, mean0, dif, dlta, 
      K = Length[T]}, type = IAType /. {opts} /. Options[IADesign]; 
      dif = Differencing /. {opts} /. Options[IADesign]; 
      dlta = Delta1 /. {opts} /. Options[IADesign]; If[K > 1 && dlta != 0, 
       Message[IADesign::notAvail]; dlta = 0; ]; If[K != Length[g], 
       Print["Error: Lengths of T and g different."]; Return[]]; 
      h = If[dlta === 0, 0, 1]; {d, ds, s} = Which[NumberQ[dif], {dif, 0, 0}, 
        VectorQ[dif], Take[dif, 3], True, {0, 0, 0}]; 
      y = Table[x = N[Which[type === Step, Table[If[i < T[[k]] + h, 0, 1], 
             {i, n + h}], type === Pulse, Table[If[i != T[[k]] + h, 0, 1], 
             {i, n + h}], type === Ramp, Table[If[i < T[[k]] + h, 0, 
              i - T[[k]] - h + 1], {i, n + h}], True, 
            Message[IADesign::invalidia, type]; Table[If[i < T[[k]] + h, 0, 
               1], {i, n + h}]]]; x = If[dlta === 0, x, 
           Drop[First /@ FoldList[Flatten[{{dlta} . #1 + #2, Drop[#1, 
                  -1]}] & , {0}, x], 1]]; x = TSDifference[x, d, ds, s]; 
         x = NestList[Drop[Prepend[#1, 0], -1] & , x, g[[k]]]; x, {k, K}]; 
      y = Flatten[y, 1]; Transpose[Prepend[y, Table[1, {Length[First[y]]}]]]];

IAFit1[(z_)?NumberVectorQ, (T_)?NumberVectorQ, (g_)?NumberVectorQ, 
     (opts___)?OptionQ] := Module[{etol, n = Length[z], \[Lambda]hat = 0., 
      zetaprev = \[Lambda]prev = Infinity, zetahat = \[Lambda]hat = 0, y, 
      res, theta = phi = 0., verb, zm, type, err, dlta, i = 0, relerr = 100., 
      K = Length[T], gkp1, om, omout, iterCount = 1, verbout, verbOpt}, 
     verbOpt = Verbose /. {opts} /. Options[IAFit1]; 
      iterCount = If[IntegerQ[verbOpt], verbOpt, 1]; 
      verb = Which[iterCount < 1, False, LogicalQ[verbOpt], verbOpt, 
        iterCount >= 1, True]; type = IAType /. {opts} /. Options[IAFit1]; 
      err = Noise /. {opts} /. Options[IAFit1]; 
      etol = ErrorTolerance /. {opts} /. Options[IAFit1]; 
      dlta = Delta1 /. {opts} /. Options[IAFit1]; 
      maxiter = MaxitIAFit1 /. {opts} /. Options[IAFit1]; 
      deltaOperator = If[dlta === 0, {}, {dlta}]; 
      res = y = z - (zm = Mean[z]); x = IADesign[n, T, g, IAType -> type, 
        Delta -> dlta]; While[maxiter > i && etol < relerr, 
       zetaprev = zetahat; \[Lambda]prev = \[Lambda]hat; 
        W = InverseCovarianceMatrixARMA[{phi}, {theta}, n]; 
        zetahat = Flatten[Inverse[Transpose[x] . W . x] . Transpose[x] . W . 
           Transpose[{y}]]; res = y - First[zetahat]; 
        zetak = Drop[zetahat, 1]; Do[gkp1 = g[[k]] + 1; 
          res -= ICalc[n, T[[k]], Take[zetak, gkp1], deltaOperator]; 
          zetak = Drop[zetak, gkp1], {k, K}]; \[Lambda]hat = 
         Which[err == "AR1", phi = AR1Fit[res, MeanZero -> True], 
          err == "MA1", theta = MA1Fit[res, MeanZero -> True], True, 0]; 
        relerr = Max[Abs /@ Flatten[{zetaprev - zetahat, \[Lambda]prev - 
               \[Lambda]hat}]/(Abs /@ Flatten[{zetahat, \[Lambda]hat}] + 
            10^(-2))]; verbout = If[Mod[i, iterCount] == 0, True, False]; 
        If[verb && verbout, Print[StyleForm[StringForm[
            "`` zetahat = ``,  \[Lambda]hat = ``  relerr = ``", i, zetahat, 
            \[Lambda]hat, relerr], FontSize -> 10, FontColor -> 
            RGBColor[0, 0, 1]]]]; If[maxiter == ++i, Message[IAFit1::maxiter, 
          i, relerr]]; If[\[Lambda]hat === 0, Break[]]; ]; 
      om = Drop[zetahat, 1]; omout = {}; 
      Do[gkp1 = g[[k]] + 1; omout = Join[omout, Take[om, gkp1]]; 
        om = Drop[om, gkp1], {k, K}]; {Mu -> zm + First[zetahat], 
       Omega -> If[K == 1, First[omout], omout], Delta -> dlta, 
       Lambda -> \[Lambda]hat, Loglikelihood -> LoglikelihoodARMA[{phi}, 
         {theta}, res, MeanZero -> True], IVType -> type, 
       IVTime -> If[K == 1, First[T], T], IVErr -> err, Residuals->res}];
 
IAFit1[(z_)?NumberVectorQ, (T_)?NumberQ, (g_)?NumberQ, (opts___)?OptionQ] := 
    Module[{etol, n = Length[z], \[Lambda]hat = 0., 
      zetaprev = \[Lambda]prev = Infinity, zetahat = \[Lambda]hat = 0, y, 
      res, theta = phi = 0., verb, zm, type, err, dlta, i = 0, 
      relerr = 100.}, verb = Verbose /. {opts} /. Options[IAFit1]; 
      type = IAType /. {opts} /. Options[IAFit1]; 
      err = Noise /. {opts} /. Options[IAFit1]; 
      etol = ErrorTolerance /. {opts} /. Options[IAFit1]; 
      dlta = Delta1 /. {opts} /. Options[IAFit1]; 
      maxiter = MaxitIAFit1 /. {opts} /. Options[IAFit1]; 
      IAFit1[z, {T}, {g}, Verbose -> verb, IAType -> type, Noise -> err, 
       ErrorTolerance -> etol, Delta1 -> dlta, MaxitIAFit1 -> maxiter]];

IAFit[(z_)?NumberVectorQ, (T_)?NumberQ, (g_)?NumberQ, (opts___)?OptionQ] := 
    Module[{n = Length[z], \[Lambda]hat = 0., 
      zetaprev = \[Lambda]prev = Infinity, zetahat = \[Lambda]hat = 0, y, y0, 
      theta = phi = 0., verb, dlta, EstDelta, dif, d}, 
     verbBrent = VerboseBrent /. {opts} /. Options[IAFit]; 
      verbIAFit1 = VerboseIAFit1 /. {opts} /. Options[IAFit]; 
      EstDelta = EstimateDelta /. {opts} /. Options[IAFit]; 
      dif = Differencing /. {opts} /. Options[IAFit]; 
      type = IAType /. {opts} /. Options[IAFit]; 
      err = Noise /. {opts} /. Options[IAFit]; If[EstDelta, 
       dlta = Brent[Loglikelihood /. IAFit1[z, T, g, Noise -> err, 
             IAType -> type, Differencing -> dif, Verbose -> verbIAFit1, 
             Delta1 -> #1] & , -0.99, 0.99, Verbose -> verbBrent]; 
        IAFit1[z, T, g, Noise -> err, IAType -> type, Differencing -> dif, 
         Verbose -> verbIAFit1, Delta1 -> dlta], IAFit1[z, T, g, 
        Noise -> err, IAType -> type, Differencing -> dif, 
        Verbose -> verbIAFit1]]];
 
TestIA[n_, T_, omega_, delta_, phi_, theta_, opts___] := 
    Module[{b = Join[omega, delta], K = Length[omega] + Length[delta], 
      \[GothicG] = Plus @@ omega/(1 - Plus @@ delta), IFull, Varb, X, xcrit1, 
      xcrit2, p1, p2, p2U, p2LZ}, sigmasquared = SigmaSquared /. {opts} /. 
        Options[TestIA]; PierceA = PierceApproximation /. {opts} /. 
        Options[TestIA]; whichTest = Test /. {opts} /. Options[TestIA]; 
      ttype = TestType /. {opts} /. Options[TestIA]; 
      iatyp = IAType /. {opts} /. Options[TestIA]; 
      IFull = InformationMatrixIA[n, T, omega, delta, phi, theta, 
         IAType -> iatyp, PierceApproximation -> PierceA]/sigmasquared; 
      InvIFull = Check[Inverse[IFull], True]; 
      Varb = If[TrueQ[InvIFull], True, Drop[Transpose[Drop[InvIFull, 1]], 
         1]]; p1 = If[TrueQ[Varb], Message[TestIA::singular]; Indeterminate, 
        X = b . Inverse[Varb] . b; 1 - CDF[ChiSquareDistribution[K], X]]; 
      Z = \[GothicG]/(seg = SEGain[omega, delta, IFull]); 
      p2U = 1 - CDF[NormalDistribution[0, 1], Z]; 
      p2L = CDF[NormalDistribution[0, 1], Z]; 
      p2 = Switch[ttype, TwoSided, 2*Min[p2U, p2L], UpperTail, p2U, 
        LowerTail, p2L, True, p2U + p2L]; Switch[whichTest, Both, 
       {WaldTest -> p1, GainTest -> p2}, GainTest, p2, Waldtest, p1, True, 
       {WaldTest -> p1, GainTest -> p2}]];

PowerIA[n_, T_, omega_, delta_, phi_, theta_, opts___] := 
 Module[{b = Join[omega, delta], K = Length[omega] + Length[delta], 
  \[GothicG] = Plus @@ omega/(1 - Plus @@ delta), IFull, Varb, X, xcrit1, 
  xcrit2, p1, p2, p2U, p2LZ}, 
 sigmasquared = SigmaSquared /. {opts} /. 
        Options[PowerIA]; 
 lev = Level /. {opts} /. Options[PowerIA]; 
 PierceA = PierceApproximation /. {opts} /. Options[PowerIA]; 
 whichTest = Test /. {opts} /. Options[PowerIA]; 
 ttype = TestType /. {opts} /. Options[PowerIA]; 
 iatyp = IAType /. {opts} /. Options[PowerIA]; 
 IFull = InformationMatrixIA[n, T, omega, delta, phi, theta, 
         IAType -> iatyp, PierceApproximation -> PierceA]/sigmasquared; 
 InvIFull = Check[Inverse[IFull], True]; 
 Varb = If[TrueQ[InvIFull], True, Drop[Transpose[Drop[InvIFull, 1]], 
         1]]; 
 p1 = If[TrueQ[Varb], Message[PowerIA::singular]; Indeterminate, 
         X = b . Inverse[Varb] . b; 
         xcrit1 = Quantile[ChiSquareDistribution[K], 1 - lev]; 
         1 - CDF[NoncentralChiSquareDistribution[K, X], xcrit1]
      ]; 
 xcrit2 = Quantile[NormalDistribution[0, 1], 1 - lev/2]; 
 Z = \[GothicG]/(seg = SEGain[omega, delta, IFull]); 
 p2U = 1 - CDF[NormalDistribution[0, 1], xcrit2 - Z]; 
 p2L = CDF[NormalDistribution[0, 1], -xcrit2 - Z]; 
 p2 = Switch[ttype, 
             TwoSided, p2U + p2L, 
             UpperTail, p2U, 
             LowerTail, p2L, 
             True, p2U + p2L]; 
 Switch[whichTest, 
        Both, {WaldTest -> p1, GainTest -> p2}, 
        Gain, p2, 
        Wald, p1, 
        True, {WaldTest -> p1, GainTest -> p2}]
];

InformationMatrixIA[n_, T_, omega_, delta_, phi_, theta_, opts___] := 
 Module[{x, g = Length[omega], h = Length[delta], Xcov, PierceA, 
         sigmasquared}, 
 sigmasquared = SigmaSquared /. {opts} /. Options[InformationMatrixIA]; 
 Xcov = X /. {opts} /. Options[InformationMatrixIA]; 
 PierceA = PierceApproximation /. {opts} /. Options[InformationMatrixIA]; 
 ityp = IAType /. {opts} /. Options[InformationMatrixIA]; 
 If[PierceA, 
  InfMatIAP[n, T, omega, delta, phi, theta, SigmaSquared -> sigmasquared], 
  x = Join[
       {Table[-1, {n}]}, 
        -Table[ICalc[n, T + i, {1}, delta, IAType -> ityp], 
        {i, 0, g - 1}]]; 
  x = If[h > 0, Join[x, 
    -Table[ICalc[n, T + i, omega, 
       -Drop[CoefficientList[ARPolynomial[delta, B]^2, B], 1], 
       IAType -> ityp], 
     {i, 1, h}]], x]; 
 x = If[Length[Xcov > 0] > 0, Join[x, Xcov], x]; 
 x . InverseCovarianceMatrixARMA[phi, theta, n] . Transpose[x]/sigmasquared]
];
 
InfMatIAP[n_, T_, omega_, delta_, phi_, theta_, opts___] := 
 Module[{sigmasquared, ityp, b, B1, b\[Mu], B\[Mu], B\[Mu]1, d, F\[Mu], 
   F1, G1, G\[Mu]}, 
 sigmasquared = SigmaSquared /. {opts} /. Options[InfMatIAP]; 
 ityp = IAType /. {opts} /. Options[InfMatIAP]; 
 b = -ICalc[n, T, Prepend[-phi, 1], 
       -Drop[
        CoefficientList[ARPolynomial[delta, B]*ARPolynomial[theta, B], B], 
        1], 
        IAType -> ityp]; 
 B1 = Toeplitz[Take[ListCorrelate[b, b, 1, 0], Length[omega]]]; 
 b\[Mu] = -TFCalc[Prepend[-phi, 1], theta, Array[1 & , {n}]]; 
 B\[Mu] = {{Plus @@ (b\[Mu]^2)}}; 
 B\[Mu]1 = {Table[ListCorrelate[b, b\[Mu], 1, 0][[j]], {j, Length[omega]}]}; 
 If[Length[delta] == 0, 
  BlockMatrix[{{B\[Mu], B\[Mu]1},{Transpose[B\[Mu]1], B1}}], 
 d = TFCalc[omega, delta, b]; 
 d\[Mu] = TFCalc[omega, delta, b\[Mu]]; 
 F\[Mu] = {Table[ListCorrelate[d, b\[Mu], 1, 0][[j]],{j, Length[omega]}]}; 
 F1 = {Table[ListCorrelate[d, b, 1, 0][[j]], {j, Length[omega]}]}; 
 G1 = Toeplitz[Take[ListCorrelate[d, d, 1, 0], Length[omega]]]; 
 G\[Mu] = {{Plus @@ (d\[Mu]^2)}}; 
 BlockMatrix[{{B\[Mu], B\[Mu]1, F\[Mu]}, {Transpose[B\[Mu]1], B1, F1}, 
           {Transpose[F\[Mu]], Transpose[F1], G1}}]/sigmasquared]
];
 
Options[InfMatIAP] = {SigmaSquared -> 1, X -> {}, 
        PierceApproximation -> False, IAType -> Step};

Options[InformationMatrixIA] = {SigmaSquared -> 1, X -> {}, 
     PierceApproximation -> False, IAType -> Step};
 
ICalc[n_, T_, omega_, delta_, (opts___)?OptionQ] := 
  Module[
   {x, g = Length[omega] - 1, h = Length[delta], type}, 
   type = IAType /. {opts} /. Options[ICalc]; 
   x = Which[
        type === Step, 
         Table[If[i < T + g, 0, 1], {i, n + g}], 
        type === Pulse, 
         Table[If[i != T + g, 0, 1],{i, n + g}], 
        type === Ramp, 
         Table[If[i < T + g, 0, i - T - g + 1], {i, n + g}], 
        True, 
         Table[If[i < T + g, 0, 1], {i, n + g}]
       ]; 
  x = Reverse[ListCorrelate[omega, Reverse[x]]]; 
  If[h == 0, x, 
   Drop[First /@ 
    FoldList[Flatten[{delta . #1 + #2, Drop[#1, -1]}] & , Table[0, {h}], x],
    1]]
  ];

TFCalc[omega_, delta_, x_] := 
 Module[{v, gp1 = Length[omega], 
    h = Length[delta], \[GothicG] = Plus @@ omega/(1 - Plus @@ delta)}, 
 v = Partition[Join[
       Table[x[[1]], {gp1 - 1}], x], gp1, 1] . Reverse[omega]; 
 If[h > 0, 
       Drop[First /@ FoldList[Flatten[{delta . #1 + #2, Drop[#1, -1]}] & , 
          Table[\[GothicG]*First[x], {h}], v], 1], 
 v]
];
 
InformationMatrixIA[n_, T_, omega_, delta_, phi_, theta_] := 
 Module[{x, g = Length[omega], h = Length[delta]}, 
 x = -Join[{Table[1, {n}]}, Table[ICalc[n, T + i, {1}, delta], {i, 0, g - 1}]]; 
 x = If[h > 0, Join[x, -Table[ICalc
        [n, T + i, omega, -Drop[CoefficientList[ARPolynomial[delta, B]^2, B], 1]], 
                         {i, 1, h}]
                    ]
     ,x]; 
 x . InverseCovarianceMatrixARMA[phi, theta, n] . Transpose[x]
];
 
VarianceGain[omega_, delta_, imatrix_] := 
 Module[{h = Length[delta], denom, g, d1, d2, d, varm}, 
  varm = Transpose[Drop[Transpose[Drop[Inverse[imatrix], 1]], 1]]; 
  denom = 1 - Plus @@ delta; 
  g = Plus @@ omega/denom; 
  d1 = 1/denom; 
  d2 = g/denom; 
  d = Table[d1, {Length[omega]}]; 
  d = If[h > 0, Join[d, Table[g/denom, {h}]], d]; 
  d . varm . d];
 
SEGain[omega_, delta_, imatrix_] := 
   Sqrt[VarianceGain[omega, delta, imatrix]];

(*=========================INPUT/OUTPUT==============================*)

ReadARMAParameters[file_] := 
Module[{pars,beta,p,q},
   pars = ReadList[file, Number];
   p = pars[[1]];
   q = pars[[2]];
   beta = Drop[pars, 2];
   Return[{Take[beta,p], Take[beta,-q]}]
   ];

ReadSARMAParameters[file_] := 
Module[{pars,beta,p,q,ps,qs,s},
   pars = ReadList[file, Number];
   p = pars[[1]];
   q = pars[[2]];
   ps = pars[[3]];
   qs = pars[[4]];
   s = pars[[5]];
   beta = Drop[pars, 5];
   Return[{Take[beta,p], Take[Drop[beta,p],q],
           Take[Drop[beta,p+q],ps],
           Take[beta,-qs], s }]
   ];
        
(*=========================Data Transformation==============================*)
TSIntegrate[(w_)?VectorQ, (d_)?IntegerQ, (ds_)?IntegerQ, (s_)?IntegerQ, 
     (zinit_)?VectorQ] := DoTSInteg[DoTSInteg[w, ds, s, 
       TSDifference[zinit, d]], d, zinit] /; d > 0 && ds > 0 && s > 0 && 
      Length[zinit] >= d + s*ds;
 
TSIntegrate[(w_)?VectorQ, (d_)?IntegerQ, (zinit_)?VectorQ] := 
    DoTSInteg[w, d, zinit] /; d > 0 && Length[zinit] >= d;
 
TSIntegrate[(w_)?VectorQ, (ds_)?IntegerQ, (s_)?IntegerQ, (zinit_)?VectorQ] := 
    DoTSInteg[w, ds, s, zinit] /; ds > 0 && s > 0 && Length[zinit] >= s*ds;
 
DoTSInteg[w_, d_, zinit_] := Module[{W = w}, 
     Do[W = Flatten[TSInteg[W, d - i + 1, Take[zinit, d - i + 1]]], {i, d}]; 
      W];
 
DoTSInteg[w_, d_, s_, zinit_] := Module[{W = w}, 
     Do[W = Flatten[TSInteg[W, d - i + 1, s, zinit]], {i, d}]; W];
 
TSInteg[w_, d_, zinit_] := FoldList[Plus, zinit, w] /; d == 1;
 
TSInteg[w_, d_, zinit_] := TSInteg[w, d - 1, 
      First[TSDifference[zinit, d - 1]]] /; d > 1;
 
TSInteg[w_, d_, s_, zinit_] := 
    Take[Flatten[FoldList[Plus, Take[zinit, s], 
        Partition[PadRight[w, s*Ceiling[Length[w]/s]], s]]], 
      Length[w] + s] /; d == 1;
 
TSInteg[w_, d_, s_, zinit_] := TSInteg[w, d - 1, s, 
                                 TSDifference[zinit,d-1,s]]; 

TSDifference[z_, d_:1 , ds_:0, s_:1] := 
TSDifference[TSDifference[z,d,1],ds,s];

TSDifference[z_, d_:1 , s_:1] := 
  With[{w = TSDifference[z, d - 1, s]}, Drop[w, s] - Drop[w, -s]];
 
TSDifference[z_, 0, s_] = z;

TSDifference[z_, d_/;d<0, s_] = Null;
 
PowerTransformation[z_, lambda_] := 
  If[Abs[lambda] < 10^(-7), Log[z], z^lambda];

TSDifference[z_, d_ /; d > 0, s_:1] := 
  With[{w = TSDifference[z, d - 1, s]}, Drop[w, s] - Drop[w, -s]];
 
TSDifference[z_, 0, s_] := z;
 
PowerTransformation[z_, lambda_] := 
  If[Abs[lambda] < 10^(-7), Log[z], z^lambda];

(*=========================ARMA Utilities==============================*)
SARMAToARMA[phi_, theta_, sphi_, stheta_, s_] := {SARToAR[phi, sphi, s], 
     SARToAR[theta, stheta, s]};
 
SARToAR[phi_, sphi_, s_] := 
    -Drop[CoefficientList[Collect[Prepend[-phi, 1] . Table[B^(i - 1), 
           {i, Length[phi] + 1}]*Prepend[-sphi, 1] . Table[B^(s (i - 1)), 
           {i, Length[sphi] + 1}], B], B], 1];
 
ARPolynomial[phi_, B_] := 
  Prepend[-phi, 1] . Table[B^j, {j, 0, Length[phi]}];

ARToMACoefficient[(phi_)?VectorQ, k_ /; k == 0] := 
  ARToMACoefficient[phi, k] = 1;
 
ARToMACoefficient[(phi_)?VectorQ, k_ /; k < 0] := 
  ARToMACoefficient[phi, k] = 0;
 
ARToMACoefficient[(phi_)?VectorQ, k_ /; k > 0] := 
  ARToMACoefficient[phi, k] = 
   Sum[phi[[i]]*ARToMACoefficient[phi, k - i], {i, 1, Length[phi]}];
 
GenerateARPolynomial[phi_] := 
  Prepend[-phi, 1] . Table[B^j, {j, 0, Length[phi]}];

GenerateARMAPolynomial[roots_List] := 
   Return[Drop[CoefficientList[Expand[-Product[1 - roots[[i]]*B, 
         {i, 1, Length[roots]}]], B], 1]];

ARPolynomial[phi_?VectorQ] := 
  Prepend[-phi, 1] . Table[B^j, {j, 0, Length[phi]}];

ARRoots[phi_] := 
         Table[
          NRoots[
                 1-Sum[phi[[i]] B^i, {i, 1, Length[phi]}] == 0, B
                ] 
          [[i,2]], {i,1,Length[phi]}
          ];

SARRoots[PHI_, S_] := Flatten[
                       Outer[Times, ARRoots[PHI], E^(2 Pi I Range[S]/S)]
                             ];

ExpandToSeasonal[phis_, s_] := 
  Flatten[(Join[Table[0, {s - 1}], {#1}] & ) /@ phis];

ARInvert[(phi_)?VectorQ, MaxLag_ /; MaxLag > 0] := 
  Module[{B}, Drop[CoefficientList[Normal[Series[1/ARPolynomial[phi, B], 
       {B, 0, MaxLag}]], B], 1]];
 
ARPolynomial[phi_, B_] := Prepend[-phi, 1] . Table[B^j, {j, 0, Length[phi]}];

ToPhiStar[(phi_)?VectorQ, d_Integer, ds_Integer, s_Integer] := 
  (-Drop[CoefficientList[Expand[ARPolynomial[phi, B]*(1 - B)^d*
          (1 - B^s)^ds], B], 1]) /;( d >= 0) && (ds >= 0) && (s > 0);
 
ImpulseCoefficientsARMA[phi_List, theta_List, MaxLag_Integer, 
   (opts___)?OptionQ] := 
  Module[{d, ds, s, gphi}, d = 
      Difference /. {opts} /. Options[ImpulseCoefficientsARMA]; 
     ds = SeasonalDifference /. {opts} /. Options[ImpulseCoefficientsARMA]; 
     s = SeasonalPeriod /. {opts} /. Options[ImpulseCoefficientsARMA]; 
     gphi = If[d == 0 && ds == 0, phi, ToPhiStar[phi, d, ds, s]]; 
     If[MaxLag > Length[theta], 
      If[Length[gphi] == 0, Join[Prepend[-theta, 1], 
        Table[0, {MaxLag - Length[theta]}]], 
       First /@ 
        FoldList[Prepend[Drop[#1, -1], #2 + gphi . #1] & , 
         Prepend[Table[0, {Length[gphi] - 1}], 1], 
         Join[-theta, Table[0, {MaxLag - Length[theta]}]]]], 
      Take[ImpulseCoefficientsARMA[gphi, theta, Length[theta] + 1], 
       MaxLag + 1]]] /; MaxLag > 0;
 
InvertedCoefficientsARMA[phi_List, theta_List, MaxLag_Integer,(opts___)?OptionQ] := 
  Module[{d, ds, s, gphi}, d = 
      Difference /. {opts} /. Options[InvertedCoefficientsARMA]; 
     ds = SeasonalDifference /. {opts} /. Options[InvertedCoefficientsARMA]; 
     s = SeasonalPeriod /. {opts} /. Options[InvertedCoefficientsARMA]; 
     gphi = If[(d == 0) && (ds == 0), phi, ToPhiStar[phi, d, ds, s]]; 
     ImpulseCoefficientsARMA[theta, gphi, MaxLag]*
      Prepend[Table[-1, {MaxLag}], 1]] /; MaxLag > 0;

Brent[func_, xmin_, xmax_, (opts___)?OptionQ] := 
    Module[{ng, etol, \[Theta], maxls, i, a, fa, b, fb, c, fc, d, fd, fx, x, 
      z = N[zz], it = 1, maxi, ver, relerr = 100}, 
     maxi = MaxIterBrent /. {opts} /. Options[Brent]; 
      etol = ErrorTolerance /. {opts} /. Options[Brent]; 
      ver = Verbose /. {opts} /. Options[Brent]; 
      ng = GridSearch /. {opts} /. Options[Brent]; 
      ls = Table[{\[Theta], func[\[Theta]]}, {\[Theta], xmin, xmax, 
         (xmax - xmin)/ng}]; lsout = ls; If[ImaginaryQ[Last[Transpose[ls]]], 
       Return[Message[Brent::imaginary]]]; 
      maxls = First[Select[ls, #1[[2]] == Max[Transpose[ls][[2]]] & ]]; 
      i = Last[Flatten[Select[Transpose[Append[Transpose[ls], 
            Range[Length[ls]]]], maxls == Take[#1, 2] & ]]]; 
      Which[i == 1, Return[xmin], i == Length[ls], Return[xmax], True, 
       {{a, fa}, {b, fb}, {c, fc}} = ls[[{i - 1, i, i + 1}]]]; 
      While[maxi > it && etol < relerr, 
       fx = func[x = b + 0.5*(((b - a)^2*(fb - fc) - (b - c)^2*(fb - fa))/
              ((b - a)*(fb - fc) - (b - c)*(fb - fa)))]; 
        {{a, fa}, {b, fb}, {c, fc}} = If[fx > fb, Which[x < b, 
           {{a, fa}, {x, fx}, {b, fb}}, x > b, {{b, fb}, {x, fx}, {c, fc}}], 
          If[b - a > 0.5*(c - a), d = b + 0.38196*(a - b); fd = func[d]; 
            If[fd > fb, {{a, fa}, {d, fd}, {b, fb}}, {{d, fd}, {b, fb}, 
              {c, fc}}], d = b + 0.38196*(c - b); fd = func[d]; 
            If[fd > fb, {{b, fb}, {d, fd}, {c, fc}}, {{a, fa}, {b, fb}, 
              {d, fd}}]]]; relerr = (c - a)/(Abs[Max[a, c]] + 0.01); 
        If[ver, Print[StyleForm[StringForm[
            "`` interval: [``,``],  relerr = ``", it, a, b, relerr], 
           FontSize -> 10, FontColor -> RGBColor[1, 0, 1]]]]; 
        If[maxi == ++it, Message[Brent::maxiter, it, relerr]]; 
        If[ImaginaryQ[relerr], Return[Message[Brent::imaginary]]]; ]; b];
 
(*=========================Utilities==============================*)
InformationMatrixToCorrelationMatrix[imatrix_] := 
 Module[{se, v}, 
  se = Sqrt /@ Transpose[v = Inverse[imatrix], {1, 1}]; 
  v/Outer[Times, se, se]
];

InformationMatrixToSE[imatrix_] := 
  Sqrt /@ Transpose[Inverse[imatrix], {1, 1}];
 
Average[x_] := With[{y=x//Flatten}, Apply[Plus, y]/Length[y]];

IndexOfMinimum[(x_)?NumberVectorQ] := 
  First[First[Select[Transpose[{Range[Length[x]], x}], 
    #1[[2]] == Min[x] & ]]];
 
RJoin[a_, b_] := 
  Which[a == {{}} && b == {{}}, {}, a == {{}} && b != {{}}, b, 
   a != {{}} && b == {{}}, a, True, 
   Transpose[Join[Transpose[a], Transpose[b]]]];
 
(*******************************************************************)

End[];
EndPackage[];
