Print[StyleForm[StringJoin[
        "Updated January 7, 2007.",
        "Loaded FitAR"
        ]
  ,FontSize->9,FontFamily->"Times", FontColor->
  RGBColor[0., 0.25,0.25]
   ]];

(* Author:  A.I. McLeod,  aimcleod@uwo.ca
	    Ying Zhang,   ying.zhang@acadiau.ca *)
(* Copyright (c), A.I. McLeod and Y. Zhang, 2005 *)
(* FitAR is a package for modelling autoregressive models *)
(* Summary:
A complete package for model building for AR processes. Topics include:
sample and partial autocorrelation functions, theoretical autocorrelation,
AR spectral estimation, exact likelihood estimation, diagnostic checks.
Comprehensive diagnostic checks include: portmanteau test, residual
autocorrelations and their consonance interval, normal probability plot,
time series plot, plot of the roots of the AR characteristic equation,
plot of the fitted theoretical autocorrelation and spectral density.
Simulation and bootstrapping of AR. Forecasting AR. Box-Cox loglikelihood.
*)

(* ORGANIZATION:
     FitAR and related functions
     diagnostic checks
     theoretical autocorrelation
     spectral density
     sample autocorrelation
     partial autocorrelation
     simulation
     utilities 
*)



Needs["Statistics`DescriptiveStatistics`"];

BeginPackage["FitAR`"]

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

(* Load need packages *)
Needs["Statistics`DescriptiveStatistics`"];
Needs["Statistics`ContinuousDistributions`"];
Needs["Statistics`DataManipulation`"];
Needs["Statistics`MultinormalDistribution`"]; 
Needs["LinearAlgebra`MatrixManipulation`"]; 

(* initializations *)
MeanZero;UseMeanMLE;Verbose;
MeanValue;
TestType; 
PlotType; 
SimultaneousBound; SquaredResiduals; FractionalDifference;
LoggedSdf; LoggedF; PlotRange; PlotLabel; Normalized; InnovationVariance;
BartlettLargeLag; Bound; PlotLabel;
Method; 
LinearPredictor;  ModelSelection;
Criterion; Table; Plot;
Silent;
ARModel; Loglikelihood; ARParameterTable; AIC;
BIC; ARCoefficients; MeanValue; 
ResidualSd; PortmanteauDiagnostic;
Pacf; Iterations; AdjustableParameters; 
Residuals; OriginalSeries;
PortmanteauStatistic; DegreesOfFreedom; Pvalue;
ARModelOrder;
ResidualVariance;
\[Lambda];
SubsetTable;MinAICModel;MinBICModel;
SubsetARQ;ExplosiveQ;AICLL;BICLL;RMSE;RSQ;ForecastTable;XMatrix;

(*FitAR and related functions*)

Options[FitARLS] = {MeanZero -> False, MaxLag -> Automatic};
 
FitARLS /: FitARLS::usage = "FitARLS[z, p] fits an autoregression or subset \
autoregression to time series z using least squares (conditional mle). If p \
is a positive integer, an AR(p) model is fit.  If p is a vector of positive \
integers, a subset \!\(AR\_\[Zeta]\)(p) model is fit using the partial \
autocorrelation parameterization.  The sample mean is used as the estimator \
of the series mean by default.  If the mean is known and subtracted from the \
data use the option MeanZero->True.  MaxLag->Automatic, the number of lags in \
the Ljung-Box portmanteau test is chosen automatically depending on the length of the \
input series. Set MaxLag->n to use exactly n lags in this test.";

Options[FitAR] = {MeanZero -> False, UseMeanMLE -> False, Verbose -> False, 
     MaxLag -> Automatic, InitialEstimate -> "Burg"};
 
FitAR /: FitAR::usage = "FitAR[z, p] fits an autoregression or subset \
autoregression to time series z. If p is a positive integer, an AR(p) model \
is fit.  If p is a vector of positive integers, a subset \!\(AR\_\[Zeta]\)(p) \
model is fit using the partial autocorrelation parameterization.  The sample \
mean is used as the estimator of the series mean by default.  If the mean is \
known and subtracted from the data use the option MeanZero->True. For \
UseMeanMLE->True, the exact mle of the mean is obtained. Set Verbose->True, \
to see a printout of the iterations when the exact mle is computed for the \
mean. MaxLag->Automatic, the number of lags in the portmanteau test is chosen \
automatically depending on the length of the input series. Set MaxLag->n to \
use exactly n lags in this test. 
InitialEstimate->\"Burg\" \
is the default but you can also set \
InitialEstimate->\"YuleWalker\",InitialEstimate->\"Zero\",InitialEstimate->\"\
Random\", or InitialEstimate->XINIT, where XINIT is a vector specifying the \
initial estimates. Note that XINIT must specify the partial autocorrelations \
(not the AR Coefficients themselves).";
 
GetARFit /: GetARFit::usage = "GetARFit[pvec,y] uses FindMinimum to minimize \
the negative of the exact loglikelihood function given SubAR lags in pvec and \
centered data vector y. It is assumed that y has mean zero. Examples to fit an AR(3), \
GetARFit[{1,2,3},y] and to fit a constrained AR(3), GetARFit[{1,3},y].  The \
default behaviour is to use the Burg estimate as an initial estimate.  To use \
something else as an initial estimate use GetARFit[pvec,y,xinit], where xinit \
specifies the vector of partial autocorrelations.";
 
Options[ChampernowneD] = {MeanZero -> False};
 
ChampernowneD /: ChampernowneD::usage = "ChampernowneD[z] returns the sample \
D matrix. It takes option MeanZero->False. If MeanZero->False, then the \
sample mean of z is subtracted and the sacvf is computed for the corrected \
series.";

MeanAR /: MeanAR::usage = "MeanAR[z,phi] computes the exact mle for the mean \
of an AR with parameter phi and data z.";
 
CoefAR /: CoefAR::usage = "CoefAR[ans], where ans is the output from \
FitARIMA, produces a table showing the fitted model parameters.";
 
SummaryAR /: SummaryAR::usage = "SummaryAR[ans], where ans is the output from \
FitAR, produces a summary of the fitted model.";

FullSummaryAR /: FullSummaryAR::usage = "FullSummaryAR[ans], where ans is the output from \
FitAR, produces a summary of the fitted model.";

Options[LoglikelihoodAR] = {MeanValue -> 0};
 
LoglikelihoodAR /: LoglikelihoodAR::usage = "LoglikelihoodAR[phi,z] evaluates \
loglikelihood function maximized over \!\(\[Sigma]\_a\%2\) in the AR(p) \
model\n,\!\(z\_t\)=\!\(\[Phi]\_1\)\!\(z\_\(t - \
2\)\)+...+\!\(\[Phi]\_p\)\!\(z\_\(t - p\)\)+\!\(a\_t\),\n where \!\(a\_t\) \
are NID(0,\!\(\[Sigma]\_a\%2\)) and \
phi={\!\(\[Phi]\_1\),...,\!\(\[Phi]\_p\)}. It takes an option MeanValue which \
has default MeanValue->0 indicating the mean is known and assumed to be zero. \
To include a mean parameter in the likelihood evaluation set \
MeanValue->\[Mu], where \[Mu] the specified value of the mean. The exact form \
of the concentrated loglikelihood is given by\n-1/2*Log[g]-n/2 Log[S/n],where \
n=Length[z], g is the covariance determinant and S is the quadratic term \
which is equal to the sum of the squared standardized prediction errors.";
 
DetAR::usage=
"DetAR[phi] computes the normalized covariance determinant of AR.";

Options[ExactAIC] = {MeanZero -> False};
 
ExactAIC /: ExactAIC::usage = "ExactAIC[ans] or ExactAIC[z,phi,k] computes \
the exact value of the AIC using the exact likelihood for a stationary AR \
model.  The argument ans is the output from FitAR or FitLSAR. Alternatively, \
the arguments z,phi,k are respectively the original times, the ar \
coefficients and k, the number of adjustable parameters. Note that the \
program substracts out the sample mean so an AR(1) has two parameters \
including the mean. If the mean is known to be zero, use the option \
MeanZero->True. This option is only used with ExactAIC[z,phi,k].";

Options[ExactBIC] = {MeanZero -> False};
 
ExactBIC /: ExactBIC::usage = "ExactBIC[ans] or ExactBIC[z,phi,k] computes \
the exact value of the AIC using the exact likelihood for a stationary AR \
model.  The argument ans is the output from FitAR or FitLSAR. Alternatively, \
the arguments z,phi,k are respectively the original times, the ar \
coefficients and k, the number of adjustable parameters. Note that the \
program substracts out the sample mean so an AR(1) has two parameters \
including the mean. If the mean is known to be zero, use the option \
MeanZero->True. This option is only used with ExactBIC[z,phi,k].";



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

BackcastResiduals /: BackcastResiduals::usage = "BackcastResiduals[z,phi,thet\
a] computes the innovations for t=1,...,n for the time series in vector z of \
length n. It is assumed that z has been mean corrected and phi and theta are \
invertible.";
 

PortmanteauTest /: PortmanteauTest::usage = "PortmanteauTest[a, M, k] \
computes the Ljung-Box Portmanteau test for residuals in \
vector a and lags in vector M given that k ARMA parameters were estimated.  \
If FitAR  or FitARLS is used to fit a model, then PortmanteauTest[ans, M], where ans is \
the output. ";
 
Options[DiagnosticPlotAR] := {PlotType -> "Racf", MaxLag -> 30};
 
DiagnosticPlotAR /: DiagnosticPlotAR::usage = "DiagnosticPlotAR[ans] produces \
one or more diagnostic plots. Both regular and subset cases are handled. The \
default is the residual autocorrelation plot with 95% joint consonance \
intervals. Other options are, PlotType->\"QQ\" for normal Q-Q plot, \
PlotType->\"Trace\" for a time series plot of the residuals, \
PlotType->\"Factors\" for plotting the reciprocal of the roots of the \
operators, PlotType->\"Tacf\" for plotting the theoretical autocorrelation \
function of the fitted model, PlotType\[Rule]\"Sdf\" for plotting the fitted \
theoretical spectral density function, PlotType\[Rule]\"All\" to produce all \
plots.";
 
Options[PlotRacf] = {MaxLag -> 30, SimultaneousBound -> True, 
     SquaredResiduals -> False, FractionalDifference -> False};
 
PlotRacf /: PlotRacf::usage = "RacfPlot[res, phi] plots the residual \
autocorrelations along with benchmark limits determined by the variance \
matrix of residual autocorrelation function for the estimated model with \
parameters given by phi. PlotRacf takes options with default values, \
Maxlag->15, SimultaneousBound->True, SquaredResiduals->False.";

Options[PlotRacfAR\[Phi]] = {MaxLag -> 30, SimultaneousBound -> True, 
     SquaredResiduals -> False};
 
PlotRacfAR\[Phi] /: PlotRacfAR\[Phi]::usage = "PlotRacfAR\[Zeta][res, \[Phi], \
k] plots the residual autocorrelations along with benchmark limits determined \
by the variance matrix of residual autocorrelation function for the estimated \
subset \!\(AR\_\[Phi]\) model with parameters \[Phi] at lags k.  Example, \
\!\(AR\_\[Phi]\)(1,3) set \[Phi]={\!\(\[Phi]\_1\),\!\(\[Phi]\_3\)) and \
k={1,3}. PlotRacf takes options with default values, Maxlag->15, \
SimultaneousBound->True, SquaredResiduals->False.";

InfMat\[Phi] /: InfMat\[Phi]::usage = "InfMat\[Phi][\[Phi],lags] returns the \
information matrix for a subset \!\(AR\_\[Phi]\) with parameters in \[Phi] \
and lags in lags.  Example, \!\(AR\_\[Phi]\)(1,3) set \
\[Phi]={\!\(\[Phi]\_1\),\!\(\[Phi]\_3\)) and lags={1,3}.";

VarianceRacfAR\[Phi] /: VarianceRacfAR\[Phi]::usage = "VarianceRacfAR\[Phi][\
\[Phi], k, m] computes the large-sample covariance matrix of the normalized \
residual autocorrelations at lags {1,2,...,m} of the subset \!\(AR\_\[Phi]\) \
model with lags in vector k and parameters in vector \[Phi].  Example, \
\!\(AR\_\[Phi]\)(1,3) set \[Phi]={\!\(\[Phi]\_1\),\!\(\[Phi]\_3\)) and \
k={1,3}. To obtain the large-sample approximation to the residual \
autocorrelation variance, use VarianceRacfAR\[Phi][\[Phi], k, m]/n, where n \
is the length of the time series.";

Options[PlotRacfAR\[Zeta]] = {MaxLag -> 30, SimultaneousBound -> True, 
     SquaredResiduals -> False};
 
PlotRacfAR\[Zeta] /: PlotRacfAR\[Zeta]::usage = "PlotRacfAR\[Zeta][res, \
\[Zeta], k] plots the residual autocorrelations along with benchmark limits \
determined by the variance matrix of residual autocorrelation function for \
the estimated subset AR\[Zeta] model with parameters \[Zeta] at lags k. PlotRacf takes \
options with default values, Maxlag->15, SimultaneousBound->True, \
SquaredResiduals->False.";

VarianceRacfAR\[Zeta] /: VarianceRacfAR\[Zeta]::usage = "VarianceRacfAR\[Zeta\
][\[Zeta], k, m] computes the large-sample covariance matrix of the \
normalized residual autocorrelations at lags {1,2,...,m} of the SubAR(k) \
model with lags in vector k and parameters in vector \[Zeta].  To obtain the \
large-sample approximation to the residual autocorrelation variance, use \
VarianceRacfAR\[Zeta][\[Zeta], k, m]/n, where n is the length of the time \
series.";
 
VarianceRacfAR /: VarianceRacfAR::usage = "VarianceRacfAR[phi, MaxLag, n] \
determines the covariance matrix for the vector of residual autocorrelations \
at lags  1, ... , MaxLag for the AR model specified by phi 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}].";
 
PlotARFactors /: PlotARFactors::usage = "PlotARFactors[phi] plots the \
reciprocal of the roots of the  polynomial specified by\n  1 - \
phi\[LeftDoubleBracket]1\[RightDoubleBracket]*B - ... - \
phi\[LeftDoubleBracket]p\[RightDoubleBracket] \!\(B\^p\),\n where \
p=Length[phi], in the complex plane along with the unit circle.";
 
TimeSeriesPlot /: TimeSeriesPlot::usage = "TimeSeriesPlot[z] \
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. Takes options Sublength, AspectRatio, PlotLabel, FrameLabel.";

Options[TimeSeriesPlot]={Sublength->200, AspectRatio->0.25, FrameLabel->None, PlotLabel->None};

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

(* THEORETICAL AUTOCORRELATION *)

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

TacvfAR /: TacvfAR::singular = "Non-invertible AR model or very close to it: \
the A Matrix is singular or ill-conditioned.";
 
TacvfAR /: TacvfAR::usage = "TacvfAR[phi, maxlag] evaluates the theoretical \
autocovariance function at lags 0,1,...,maxlag for the AR model specified \
with autoregressive component, phi.";
 

AcfPlot /: 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.";

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.";

InformationMatrixAR::usage=
"InformationMatrixAR[phi] computes the large-sample information per observation
in an autoregression with parameters phi.";

InformationMatrixZeta::usage = "InformationMatrixZeta[\[Zeta]] or \
InformationMatrixZeta[\[Zeta],lags] computes the Fisher information matrix for AR or \
subset AR with \[Zeta] parameter";

InfMat\[Zeta]::usage = "InfMat\[Zeta][\[Zeta]] or \
InfMat\[Zeta][\[Zeta],lags] computes the Fisher information matrix for AR or \
subset AR with \[Zeta] parameter";
 
Jacobian\[Phi]\[Zeta] /: Jacobian\[Phi]\[Zeta]::usage = "Compute the Jacobian \
of the transformation \
\!\(TraditionalForm\`\[Phi]\[LongLeftRightArrow]\[Zeta], \\\\n \
\(TraditionalForm\`\[PartialD]\[Phi]\)/\(TraditionalForm\`\[PartialD]\[Zeta]\
\)\)";
  
\[GothicCapitalJ] /: \[GothicCapitalJ]::usage = "\[GothicCapitalJ][\[Zeta],k] \
k-th term in Jacobian computation\!\(TraditionalForm\`\\\\n \
\(TraditionalForm\`\[PartialD]\[Phi]\)/\(TraditionalForm\`\[PartialD]\[Zeta]\
\)\)"; 
 
(* SPECTRAL DENSITY FUNCTION *)

Options[SdfARPlot] = {LoggedSdf -> True, LoggedF -> False, PlotRange -> All, 
     PlotLabel -> None, Normalized -> False, InnovationVariance -> 1, RadialFrequency->True };
 
SdfARPlot /: SdfARPlot::usage = "SdfARPlot[phi] plots the spectral density \
function of AR. For options see Options[SdfARPlot].";
 
Options[SdfAR] = {InnovationVariance -> 1, Normalized -> False};
 
SdfAR /: SdfAR::usage = "SdfAR[phi] evaluates the spectral density function \
of an  autoregressive model,\n \
\!\(z\_t\)=\!\(a\_t\)+\!\(\[Phi]\_1\)\!\(z\_\(t - \
1\)\)+...+\!\(\[Phi]\_p\)\!\(z\_\(t - p\)\)\n at frequencies \
0,1/1024,...,511/1024 using the FFT, \
phi={\!\(\[Phi]\_1\),...,\!\(\[Phi]\_p\)}. SdfAR[phi,\[Lambda]] evaluates the \
sdf at frequency \[Lambda].";
 
(* SAMPLE AUTOCORRELATION *)

Options[SacfPlot] = {MeanZero -> False, MaxLag -> 15, BartlettLargeLag -> 0, 
     Bound -> "Individual", PlotLabel -> None};
 
SacfPlot /: 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[Sacf] = {MaxLag -> 1, MeanZero -> False};
 
Sacf /: 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.";
 
Options[Sacvf] = {MaxLag -> 1, MeanZero -> False};
 
Sacvf /: 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.";
 
(* PARTIAL AUTOCORRELATION *)

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

SPacfPlot::usage=
"SPacfPlot[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[PacfPlot] = {MeanZero -> False, MaxLag -> 15, NumberOfPanels -> 1, 
     Largest -> 0, Method -> Automatic}
 
PacfPlot /: PacfPlot::usage = "PacfPlot[z] plots the 95% confidence intervals \
for the partial autocorrelations of time series z. The partial \
autocorrelations are computed using the Burg algorithm.Options and their \
defaults are MeanZero->False, MaxLag->15, NumberOfPanels->1, Largest->0. If \
the mean of the input series has been subtracted off the input then set \
MeanZero\[Rule]True.  MaxLag controls the number of lags plotted.  When \
MaxLag is larger than 30 to 50, it is convenient to plot several panels.  The \
number of panels is controlled by NumberOfPanels. By default with \
Method->Automatic, if MaxLag is less than or equal to 30 exact mle is used \
and otherwise the Burg algorithm is used.  Otherwise, if Method->\"MLE\", \
then exact mle will be used; if Method->\"Burg\" will be used; if \
Method->\"YW\", Yule-Walker is used."
 
 
Options[PacfBurg] = {MeanZero -> False, LinearPredictor -> False, 
     ModelSelection -> False};
 
PacfBurg /: 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.";
 
Options[ICFromPacf] = {Criterion -> "AIC", Table -> False, Plot -> False};
 
ICFromPacf /: 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->k a vector of length k \
is returned contain the best k IC values.";
 
PacfToAR /: PacfToAR::usage = "PacfToAR[pi] converts the pacf in pi to \
coefficients in an AR(p), p=Length[pi]";

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

Options[PacfDL] = {LinearPredictor -> False};
 
PacfDL /: 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[AICSubAR] = {MeanZero -> False, MaxLag -> Automatic, Largest -> 5, 
     Method -> Automatic};
 
AICSubAR /: AICSubAR::usage = "AICSubAR[z] optimal AIC and BIC subset AR \
models are found based.  The pacf is generated for lags 1,...,MaxLag and the \
K values with largest absolute magnitude are selected.  The option \
MaxLag->Automatic, means that Min[n/4, 200] lags will be used by default.  By \
default Largest->5 and this is used for K.  A table is computed showing the \
AIC and BIC values.  Plots are produced.  The automatic has the form \
{graphics--, --graphics, SubsetTable->t, MinAICModel->a, MinBICModel->b}.  To \
improve the resolution of the plots the AIC and BIC values corresponding to \
the white noise model, AR(0), is not shown since by definition it is 0. As \
with PacfPlot, Method->Automatic means exact mle is used when MaxLag<=30 and \
otherwise Burg.";

(* SIMULATION *)

SimulateGaussianAR /: SimulateGaussianAR::usage = "SimulateGaussianAR[phi, n] \
simulates n consecutive observations from a stationary-causal Gaussian AR(p) \
process,\n \!\(z\_t\)=\!\(a\_t\)+\!\(\[Phi]\_1\)\!\(z\_\(t - \
1\)\)+...+\!\(\[Phi]\_p\)\!\(z\_\(t - p\)\)\n where \!\(a\_t\) are NID(0,1) \
and phi={\!\(\[Phi]\_1\),...,\!\(\[Phi]\_p\)}. For AR(1), simply use, \
phi=\!\(\[Phi]\_1\). For a time series fit with ans=ARFit, \
SimulateGaussianAR[ans] generates a simulation of the fitted model.";

BootAR /: BootAR::usage = 
     "BootAR[obj] parametric bootstrap of model fitted by FitAR or FitARLS.";
 

(* FORECAST *)
 
ForecastOneStepTestData /: ForecastOneStepTestData::usage = "ForecastOneStepT\
estData[ans,yTest] produces one-step forecasts from model contained in ans, \
the output from FitAR, and compares with out-of-sample data in yTest.";



(* UTILITIES *)

Qd::usage="Qd[phi,D] fast evaluation of phi.D.phi.";

IndexOfMinimum /: 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.";

NumberVectorQ::usage=
  "NumberVectorQ[x] returns True if x is a numeric vector.";

InvertibleQ /: 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.";

ToeplitzQ::usage=
  "ToeplitzQ[x] tests if x is a Toeplitz matrix.";

Options[ReadData] = {Silent -> False};
 
ReadData /: ReadData::NotFound = "``";
 
ReadData /: ReadData::usage = "ReadData[file] file is a character string \
which specifies a file. ReadData[file] inputs a  data file. Initial lines in \
the file which begin with the character # are assumed to contain comments and \
documentation. The remainder of the file contains the data values \
\!\(z\_1\),\!\(z\_2\), ... which are read by rows.  The output is a vector \
containing all the data. There is one option which has default value, \
Silent->False.  When Silent->False, the title is displayed as well as the \
number of data values input. When Silent->True, the title is not displayed \
nor the number of data values input.  If the file contains p columns and n \
rows, you can obtain the n-by-p matrix by using Partition.  For example, \
Partition[ReadData[file],p] returns an n-by-p matrix assuming that there are \
n*p data values present in the file.";


ARInvert::usage = "ARInvert[phi, MaxLag] computes the coefficients \
psi[1],...,psi[M] in the expansion 1/phi(B)=1+psi[1]B+...+psi[M]B^M, \
where M=MaxLag.";

TakeMatrixColumns /: TakeMatrixColumns::usage = "TakeMatrixColumns[J, i] \
returns matrix obtained by selecting the columns in vector i from matrix J.";

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

BICFromPacf /: BICFromPacf::usage = "BICFromPacf[pacf,n,k] BIC as a function \
of pacf, n=Length of series, k=# parameters";
			    
AICFromPacf /: AICFromPacf::usage = "AICFromPacf[pacf,n,k] AIC as a function \
of pacf, n=Length of series, k=# parameters";


(********************************************************************) 
Begin["`Private`"]
(********************************************************************) 

(* Least squares and subset least squares *)

MeanCorrect = Compile[{{X, _Real, 2}}, 
    With[{n = Dimensions[X] // First}, X - Table[Total /@ (Transpose[X])/n, {n}]]];

FitARLS[(y_)?NumberVectorQ, (pvec_)?NumberVectorQ, (opts___)?OptionQ] := 
    Module[{n = Length[y], pmax = Max[pvec], indmean, ym, nef, y0, X, LL, k, 
      res, rv, ise, b, resSd, se\[Beta], sdmn, se, tb, indSubsetQ, ml, ty}, 
     indmean = MeanZero /. {opts} /. Options[FitARLS]; 
      ml = MaxLag /. {opts} /. Options[FitARLS]; 
      ml = If[ml === Automatic, Which[n <= 100, 15, n <= 200, 20, n <= 300, 
         25, n <= 500, 40, True, 50], ml]; 
      ty = TestType /. {opts} /. Options[FitARLS]; 
      ym = If[indmean, 0, Mean[y]]; nef = n - pmax; 
      kAR = Length[Select[pvec, #1 != 0 & ]]; k = If[indmean, kAR, kAR + 1]; 
      y0 = Drop[N[y], pmax]; y0 = y0 - If[indmean, 0, Mean[y0]]; 
      X = Transpose[Table[Drop[Drop[N[y], pmax - pvec[[i]]], -pvec[[i]]], 
         {i, Length[pvec]}]]; 
      X = If[indmean, X, MeanCorrect[X]]; 
      b = LinearSolve[Transpose[X] . X, Transpose[X] . y0]; 
      \[Phi]s = Array[0 & , pmax]; \[Phi]s[[pvec]] = b; 
      ise =  !InvertibleQ[b]; pf = If[ise, "Not Available", 
        PacfDL[With[{r = TacvfAR[b, pmax]}, Drop[r, 1]/r[[1]]]]]; 
      res = y0 - X . b; rv = Total[res^2]/Length[res]; 
      LL = LoglikelihoodAR[\[Phi]s, y - Mean[y]]; resSd = Sqrt[rv]; 
      Inf\[Beta] = Inverse[SiddiquiMatrix[b, pmax]]; 
      se\[Beta] = Sqrt[Transpose[Inverse[Inf\[Beta]]/n, {1, 1}][[pvec]]]; 
      sdmn = resSd*Sqrt[1/n]*Sqrt[1/Plus @@ Flatten[{1, -b}]^2]; 
      se = Append[se\[Beta], sdmn]; b = Join[b, {ym, resSd}]; 
      indSubsetQ = kAR == Length[pvec]; 
      trows = Join[If[pmax > 0, Table[Format[Subscript["\[Phi]", 
            ToString[pvec[[i]]]]], {i, Length[pvec]}], {}], 
        {"\[Mu]", "\!\(\[Sigma]\_a\)"}]; 
      tb = TableForm[Transpose[{b, Append[se, " "], Append[Drop[b, -1]/se, 
           " "]}], TableHeadings -> {trows, {"Parameter", "SE", "Z"}}]; 
      pt = PortmanteauTest[res, ml, kAR]; 
      {ARModel -> pvec, Loglikelihood -> LL, ARParameterTable -> tb, 
       AIC -> AICFromPacf[pf, n, k], AICLL -> ExactAIC[y, \[Phi]s, k], 
       BIC -> BICFromPacf[pf, n, k], BICLL -> ExactBIC[y, \[Phi]s, k], 
       ARCoefficients -> \[Phi]s, MeanValue -> ym, ResidualSd -> resSd, 
       PortmanteauDiagnostic -> pt, Pacf -> ARToPacf[\[Phi]s], AdjustableParameters -> k, 
       ExplosiveQ -> ise, SubsetARQ -> indSubsetQ, MeanZero -> indmean, 
       Residuals -> res, OriginalSeries -> y}];

FitAR[(y_)?NumberVectorQ, (p_)?(NumberVectorQ[#1] || NumberQ[#1] & ), 
     (opts___)?OptionQ] := Module[{indmean, indReport, indUseMLE, indVerb, 
      ml, mu, ycor, ycorAdj, pvec, pmax, pnot, \[DoubleStruckCapitalD], n, 
      zinit, zta, ztaHat, zeps, vars, ans, phis, muHatPrev, muHat = 0, beta, 
      trows, tb, se\[Phi], In\[Phi], cor\[Phi], sdmn, LL, k, ai, bi, res, 
      resSd, pt, MaxIter = 100, ErrMu = 1, ErrLL = 1, ETOL = 10^(-6), 
      iter = 0, LLPrev, indSubsetQ, kAR, Cov\[Beta] }, 
     indmean = MeanZero /. {opts} /. Options[FitAR]; 
      indUseMLE = UseMeanMLE /. {opts} /. Options[FitAR]; 
      indVerb = Verbose /. {opts} /. Options[FitAR]; 
      ml = MaxLag /. {opts} /. Options[FitAR]; 
      ty = TestType /. {opts} /. Options[FitAR]; n = Length[y]; 
      initq = InitialEstimate /. {opts} /. Options[FitAR]; n = Length[y]; 
      ml = If[ml === Automatic, Which[n <= 100, 15, n <= 200, 20, n <= 300, 
         25, n <= 500, 40,  n<=1000, 100, n<=2000, 200, True, 300], ml]; 
      mu = If[indmean, 0, Mean[y]]; 
      pvec = If[VectorQ[p], Sort[p], Range[p]]; pmax = Max[pvec]; 
      indSubsetQ = pmax != Length[pvec]; ycor = y - mu; 
      xinit = If[StringQ[initq], Which[StringTake[initq, {1}] === "Y", 
         PacfDL[Sacvf[y, MaxLag -> pmax]], StringTake[initq, {1}] === "Z", 
         Array[0 & , pmax], StringTake[initq, {1}] === "R", 
         Table[Random[Real, {-0.25, 0.25}], {pmax}], True, 
         PacfBurg[y, pmax]], PacfBurg[y, pmax]]; 
      xinit = If[ !StringQ[initq] && NumberVectorQ[initq] && 
         Length[initq] == pmax && Max[Abs[initq]] < 1, initq, xinit]; 
      xinit = Table[xinit[[pvec[[i]]]], {i, Length[pvec]}]; 
      ans = GetARFit[pvec, ycor, xinit]; ycorAdj = ycor; muHatPrev = mu + 1.; 
      LLPrev = First[ans]; If[indUseMLE, ztaHat = Last[ans]; 
        phis = PacfToAR[ztaHat]; muHatPrev = 0; 
        muHat = MeanAR[ycorAdj, phis]; ycorAdj = ycor - muHat; 
        While[++iter < MaxIter && Max[ErrMu, ErrLL] > ETOL, 
         ans = GetARFit[pvec, ycorAdj]; ztaHat = Last[ans]; 
          phis = PacfToAR[ztaHat]; muHatPrev = muHat; 
          muHat = MeanAR[ycor, phis]; ycorAdj = ycor - muHat; 
          ErrMu = Abs[muHat - muHatPrev]/(Abs[muHat] + 0.0001); 
          ErrLL = Abs[First[ans] - LLPrev]/Abs[LLPrev]; LLPrev = First[ans]; 
          If[indVerb, Print[StyleForm[StringForm[
              "Iteration: ``,  muhat correction = `` , ErrMu = ``, ErrLL = \
``", iter, muHat, ErrMu, ErrLL], FontFamily -> Helvetica, FontColor -> 
              RGBColor[1, 0.5, 1], FontSize -> 12, FontWeight -> 
              "Bold"]]]; ]]; ztaHat = Last[ans]; phis = PacfToAR[ztaHat]; 
      If[ !indReport, Return[phis]]; mu = mu + muHat; 
      resSd = Sqrt[Qd[Prepend[PacfToAR[ztaHat], -1.], ChampernowneD[y - mu, 
           pmax, MeanZero -> True]]/n]; 
Cov\[Beta] = 
       If[indSubsetQ, Inverse[InformationMatrixZeta[ztaHat[[pvec]], pvec]], 
         SiddiquiMatrix[phis, pmax]];
 se\[Beta] = 
       Sqrt[Transpose[Cov\[Beta]/n, {1, 1}]]; 
      sdmn = resSd*Sqrt[1/n]*Sqrt[1/Plus @@ Flatten[{1, -phis}]^2]; 
      se = Append[se\[Beta], sdmn]; beta = If[indSubsetQ, 
        Flatten[{ztaHat[[pvec]], mu, resSd}], Flatten[{phis, mu, resSd}]]; 
      trows = Join[If[pmax > 0, Table[Format[Subscript[If[indSubsetQ, 
             "\[Zeta]", "\[Phi]"], ToString[pvec[[i]]]]], {i, Length[pvec]}], 
         {}], {"\[Mu]", "\!\(\[Sigma]\_a\)"}]; 
      tb = TableForm[Transpose[{beta, Append[se, " "], 
          Append[Drop[beta, -1]/se, " "]}], TableHeadings -> 
         {trows, {"Parameter", "SE", "Z"}}]; 
      LL = LoglikelihoodAR[phis, y - mu]; k = Length[pvec]; 
      kAR = Length[Select[pvec, #1 != 0 & ]]; k = If[indmean, kAR, kAR + 1]; 
      ai = -2*LL + 2*k; bi = -2*LL + Log[n]*k; 
      res = BackcastResiduals[y - mu, phis, {}];
      pt = PortmanteauTest[res, ml, kAR]; 
      {ARModel -> If[indSubsetQ, pvec, p], Loglikelihood -> LL, 
       ARParameterTable -> tb, AIC -> AICFromPacf[ztaHat, n, k], AICLL -> ai, 
       BIC -> BICFromPacf[ztaHat, n, k], BICLL -> bi, ARCoefficients -> phis, 
       MeanValue -> mu, ResidualSd -> resSd, PortmanteauDiagnostic -> pt, 
       Pacf -> ztaHat, Iterations -> iter, AdjustableParameters -> k, 
       MeanZero -> indmean, SubsetARQ -> indSubsetQ, Residuals -> res, 
       OriginalSeries -> y}];
 

GetARFit[(pvec_)?NumberVectorQ, (y_)?NumberVectorQ] := 
    Module[{zinit, pmax = Max[pvec]}, zinit = PacfBurg[y, pmax]; 
      GetARFit[pvec, y, zinit]];
 
GetARFit[(pvec_)?NumberVectorQ, (y_)?NumberVectorQ, 
     (xinit_)?NumberVectorQ] := Module[{\[DoubleStruckCapitalD], zinit, zeps, 
      vars, n = Length[y], peff = Length[pvec], pmax = Max[pvec]}, 
     zta = Table[If[MemberQ[pvec, i], Subscript[\[ScriptZ], i], 0], 
        {i, pmax}]; \[DoubleStruckCapitalD] = ChampernowneD[y, pmax, 
        MeanZero -> True]; zinit = xinit; zeps = Min[1 - Abs[zinit]]*0.5; 
      vars = Table[{Subscript[\[ScriptZ], pvec[[i]]], -zeps + zinit[[i]], 
         zeps + zinit[[i]], -0.99999, 0.99999}, {i, peff}]; 
      ans = Check[FindMinimum[Hold[0.5*Log[DtaFUN[zta]] + 
           0.5*n*Log[Qd[Prepend[PacfToAR[zta], -1.], \[DoubleStruckCapitalD]]/
              n]], Evaluate[Sequence @@ vars], Compiled -> False], 
        {Infinity, xinit}]; 
      If[First[ans] < Infinity, {First[ans], zta /. Last[ans]}, ans]];
 
ChampernowneD[z_, p_, (opts___)?OptionQ] := 
    Module[{n = Length[z], ifmean0, X, x0, x, c, \[GothicCapitalE], A}, 
     If[Plus @@ Abs[z] == 0, Print["constant series!"]; Return[{}]]; 
      ifmean0 = MeanZero /. {opts} /. Options[ChampernowneD]; 
      X = {x0 = x = z - If[ifmean0, 0, Plus @@ z/Length[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][[j,i]] = \[GothicCapitalE][[i,j]], {j, p + 1}, 
       {i, j}]; A - \[GothicCapitalE]];
 
MeanAR[z_, phi_] := Module[{g1}, g1 = get1G[phi, Length[z]]; 
      Plus @@ (g1*z)/Plus @@ g1];
 
get1G[phi_, n_] := Module[{p = Length[phi]}, 
     x0 = Plus @@ Flatten[{1, -phi}]^2; 
      x = x0 - Apply[Plus, FromSymmetricStorageUpper[
          Flatten[Table[ListCorrelate[Drop[phi, i], Flatten[{Drop[phi, i], 
               Array[0 & , p - i - 1]}]], {i, 0, p - 1}]]], {1}] - 
        get\[Kappa][phi]; Flatten[{x, Table[x0, {i, n - 2*p}], Reverse[x]}]];
 
FromSymmetricStorageUpper[x_] := Module[{n}, 
     n = Floor[(-1 + Sqrt[1 + 8*Length[x]])/2]; 
      Table[If[i <= j, x[[n*(i - 1) + j - i*((i - 1)/2)]], 
        x[[n*(j - 1) + i - j*((j - 1)/2)]]], {j, n}, {i, n}]];
 
get\[Kappa][phi_] := Reverse[Drop[FoldList[Plus, 0, 
       Reverse[Drop[TacvfMA[phi, Length[phi]], 1]]], 1]];
 
TacvfMA[(theta_List)?NumberVectorQ, MaxLag_Integer] := 
    Module[{q, b, c, rp1}, If[ !And @@ InvertibleQ /@ {theta, {}}, 
       Message[TacvfMA::singular]; Return[{}]]; q = Length[theta]; 
      rp1 = 1 + q; If[q == 0, Prepend[Table[0., {MaxLag}], 1], 
       c = getc[1, q, {0.}, theta]; b = getb[theta, c, rp1, rp1]; 
        g = getg[{0.}, b, 1, rp1]; getr[{0.}, g, 1, rp1, MaxLag]]];
  
CoefAR[ans_] := Module[{}, {ARParameterTable /. ans}//First];
 
Options[LoglikelihoodAR] = {MeanValue -> 0};
 
SummaryAR[ans_] := Module[{}, {Loglikelihood -> (Loglikelihood /. ans), 
      AIC -> (AIC /. ans), BIC -> (BIC /. ans), ARModelOrder -> 
       (ARModel /. ans), PortmanteauDiagnostic -> (PortmanteauDiagnostic /. 
        ans)}];

FullSummaryAR[ans_]:=Drop[Drop[ans,-2],{3}];
 
LoglikelihoodAR[(phi_)?VectorQ, (z_)?VectorQ, (opts___)?OptionQ] := 
    Module[{p = Length[phi], phis, n = Length[z], \[Mu]}, 
     phis = Flatten[{1, -phi}]; \[Mu] = MeanValue /. {opts} /. 
        Options[LoglikelihoodAR]; FullSimplify[-Log[DetAR[phi]]/2 - 
        n*(Log[phis . ChampernowneD[z - \[Mu], p, MeanZero -> True] . phis/n]/
          2)]];
 
DetAR[(phi_)?VectorQ] := With[{z = ARToPacf[phi]}, 
     Times @@ ((1 - z^2)^(-Range[Length[z]]))];

ExactAIC[(ans_)?ListQ] := Module[{phi, n, k, z, mQ, mz}, 
     phi = ARCoefficients /. ans; n = Length[OriginalSeries /. ans]; 
      k = AdjustableParameters /. ans; z = OriginalSeries /. ans; 
      mQ = MeanZero /. ans; mz = If[mQ, 0, Mean[z]]; 
      -2*LoglikelihoodAR[phi, z - mz] + k*2];
 
ExactAIC[(z_)?NumberVectorQ, (phi_)?NumberVectorQ, (k_)?NumberQ, 
     (opts___)?OptionQ] := Module[{n, mQ, mz}, n = Length[z]; 
      mQ = MeanZero /. {opts} /. Options[ExactAIC]; mz = If[mQ, 0, Mean[z]]; 
      -2*LoglikelihoodAR[phi, z - mz] + k*2];

ExactBIC[(ans_)?ListQ] := Module[{phi, n, k, z, mQ, mz}, 
     phi = ARCoefficients /. ans; n = Length[OriginalSeries /. ans]; 
      k = AdjustableParameters /. ans; z = OriginalSeries /. ans; 
      mQ = MeanZero /. ans; mz = If[mQ, 0, Mean[z]]; 
      -2*LoglikelihoodAR[phi, z - mz] + k*Log[n]];
 
ExactBIC[(z_)?NumberVectorQ, (phi_)?NumberVectorQ, (k_)?NumberQ, 
     (opts___)?OptionQ] := Module[{n, mQ, mz}, n = Length[z]; 
      mQ = MeanZero /. {opts} /. Options[ExactBIC]; mz = If[mQ, 0, Mean[z]]; 
      -2*LoglikelihoodAR[phi, z - mz] + k*Log[n]];
 
(* DIAGNOSTIC CHECKS *)

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]]}]];

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}}]}}]];

 
BackcastResiduals[(z_)?NumberVectorQ, (phi_)?NumberVectorQ, 
     (theta_)?NumberVectorQ] := Module[{e, thstar, zb, a, u, v, Q = 100, p, 
      q, n, \[Phi], \[Theta]}, \[Phi] = If[phi === {}, {0.}, phi]; 
      \[Theta] = If[theta === {}, {0.}, theta]; p = Length[\[Phi]]; 
      q = Length[\[Theta]]; n = Length[z]; 
      e = Join[Table[0, {q - 1}], With[{\[Beta] = Join[-\[Phi], \[Theta]], 
          y = Reverse[z], p = Length[\[Phi]], q = Length[\[Theta]]}, 
         (#1[[p + 1]] & ) /@ FoldList[Flatten[{#2, Drop[#1, {p, p + q}], 
              \[Beta] . #1 + #2, Take[#1, {p + 1, p + q - 1}]}] & , 
           Join[Reverse[Take[y, p]], Table[0, {q}]], Drop[y, p]]]]; 
      thstar = \[Theta]; u = -With[{q = Length[\[Theta]]}, 
         Prepend[Table[thstar = Append[Drop[thstar, 1], 0], {q - 1}], 
           \[Theta]] . Reverse[Take[e, -q]]]; u = Join[u, Table[0, {Q - q}]]; 
      zb = Rest[First /@ FoldList[Flatten[{\[Phi] . #1 + #2, 
             Drop[#1, -1]}] & , Take[z, p], u]]; 
      v = Partition[Join[Table[0, {p}], Reverse[zb], z], p + 1, 1] . 
        Reverse[Prepend[-\[Phi], 1]]; 
      a = First /@ FoldList[Flatten[{\[Theta] . #1 + #2, Drop[#1, -1]}] & , 
         Table[0, {Length[\[Theta]]}], v]; Take[a, -n]];

PortmanteauTest[(a_)?NumberVectorQ, (M_)?IntegerQ /; M > 0,(k_)?IntegerQ /; k >= 0] := 
  Module[{n = Length[a], r, Q, df, pv, MK, M0, sr}, 
      M0 = M;
      r = Sacf[a, MaxLag -> M0, MeanZero -> True]; 
      Q = n*(n + 2)*Plus @@ (r*(r/(n - Range[M0])));
      df = M0 - k; 
      pv = 1 - CDF[ChiSquareDistribution[df], Q];
      {PortmanteauStatistic -> Q, MaxLag -> M0, Pvalue -> pv}
      ];
 
PortmanteauTest[(ans_)?ListQ, (M_)?IntegerQ /; M > 0] := 
 Module[{model, a, k}, 
 a = Residuals /. ans; 
 model = ARModel /. ans;
 k = If[IntegerQ[model], model, Length[model]]; 
 PortmanteauTest[a, M, k]
];
 
 
DiagnosticPlotAR[(ans_)?ListQ, (opts___)?OptionQ] := 
 Module[{ty, mxl, g\[Phi], g, sQ}, 
   ty = PlotType /. {opts} /. Options[DiagnosticPlotAR]; 
   g\[Phi] = ARCoefficients /. ans; 
   mxl = MaxLag /. {opts} /.Options[DiagnosticPlotAR]; 
   mxl = If[mxl < Length[g\[Phi]],2*Length[g\[Phi]], mxl]; 
   sQ = SubsetARQ /. ans;
   Which[
      ty == "Racf", If[sQ /. ans, 
      lags = ARModel /. ans;
      PlotRacfAR\[Zeta][Residuals /. ans, (Pacf /. ans)[[lags]], lags, MaxLag -> mxl], 
      PlotRacf[Residuals /. ans,g\[Phi], MaxLag -> mxl]], 
      ty == "Factors", PlotARFactors[g\[Phi]], 
       ty == "Trace", TimeSeriesPlot[Residuals /. ans], ty == "QQ", 
       NPlot[Residuals /. ans], ty == "Tacf", g = TacvfAR[g\[Phi], mxl]; 
        AcfPlot[Drop[g, 1]/First[g]], ty == "Sdf", SdfARPlot[g\[Phi], 
        InnovationVariance -> ResidualSd /. ans], True, 
       If[sQ /. ans, lags = ARModel /. ans; PlotRacfAR\[Zeta][
           Residuals /. ans, (Pacf /. ans)[[lags]], lags], 
         PlotRacf[Residuals /. ans, g\[Phi]]]; 
	 If[Length[g\[Phi]]<=20,PlotARFactors[g\[Phi]]]; 
        TimeSeriesPlot[Residuals /. ans]; NPlot[Residuals /. ans]; 
        g = TacvfAR[g\[Phi], mxl]; AcfPlot[Drop[g, 1]/First[g]]; 
        SdfARPlot[g\[Phi], InnovationVariance -> ResidualSd /. ans]]; ];
 
PlotRacf[res_, phi_, (opts___)?OptionQ] := 
    Module[{maxlag, acfrange, frameticks, r, blim, sbound, squared, indfdQ, 
      fd}, maxlag = MaxLag /. {opts} /. Options[PlotRacf]; 
      sbound = SimultaneousBound /. {opts} /. Options[PlotRacf]; 
      indNullModel = If[Length[phi] == 0, True, False]; 
      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 || indNullModel, Table[1/Sqrt[Length[res]], {maxlag}], 
        Sqrt[Transpose[VarianceRacfAR[phi, 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]}; 
      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}], Chop[FrameTicks -> frameticks]]]];

VarianceRacfAR[phi_, MaxLag_, n_] := Module[{X = {}, x, p = Length[phi]}, 
     x = If[p > 0, -Prepend[ARInvert[phi, MaxLag - 1], 1], {}]; 
      Do[AppendTo[X, Join[Table[0, {i}], Drop[x, -i]]], {i, 0, p - 1}]; 
      X = Simplify[(IdentityMatrix[MaxLag] - Transpose[X] . 
           Inverse[InformationMatrixAR[phi]] . X)/n]];


TimeSeriesPlot[(z_)?NumberVectorQ,  opts___?OptionQ] := 
  Module[{z2, numpanels, itake, g, sl, as, fl, pl}, 
  sl = Sublength/.{opts}/.Options[TimeSeriesPlot];
  as = AspectRatio/.{opts}/.Options[TimeSeriesPlot];
  fl = FrameLabel/.{opts}/.Options[TimeSeriesPlot];
  pl = PlotLabel/.{opts}/.Options[TimeSeriesPlot];
  numpanels = If[sl >= Length[z], 1, 
     z2 = Partition[Join[Transpose[{Range[Length[z]], z}], 
                    Table[0, {sl}]], sl, sl - 1]; 
     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, FrameLabel->fl, 
             PlotRange -> {All, {Min[z], Max[z]}}, Frame -> True, 
             Axes -> False, AspectRatio -> as, DisplayFunction -> 
              Identity] & ) /@ z2; 
   Show[GraphicsArray[Partition[g, 1],PlotLabel->pl ]], 
   ListPlot[z, PlotJoined -> True, PlotRange -> {All, {Min[z], Max[z]}}, 
         Frame -> True, Axes -> False, AspectRatio -> as,  FrameLabel->fl,PlotLabel->pl] 
           ]
];

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]]}}]];
 
PlotARFactors[phi_] := Module[{framelab}, 
     framelab = FontForm[ToString[(NumberForm[#1, 2] & ) /@ phi], 
        {"Helvetica-Bold", 10}]; 
      Show[Graphics[{{RGBColor[1, 1, 0], Disk[{0, 0}, 1, {-Pi, Pi}]}, 
         {RGBColor[0, 0, 1], PointSize[0.02], Point /@ 
           (Transpose[{Re[#1], Im[#1]}] & )[1/ARRoots[phi]]}}, 
        AspectRatio -> Automatic, Frame -> True, Axes -> True, 
        AxesOrigin -> {0, 0}, AxesLabel -> {"Real", "Imaginary"}, 
        GridLines -> {Range[-1, 1, 0.1], Range[-1, 1, 0.1]}, 
        FrameLabel -> {framelab, ""}]]];

PlotRacfAR\[Zeta][res_, \[Zeta]_, lags_, (opts___)?OptionQ] := 
    Module[{maxlag, acfrange, frameticks, r, blim, sbound, squared, n}, 
     n = Length[res]; maxlag = MaxLag /. {opts} /. 
        Options[PlotRacfAR\[Zeta]]; maxlag = If[maxlag <= Max[lags], 
        2*Max[lags], maxlag]; sbound = SimultaneousBound /. {opts} /. 
        Options[PlotRacfAR\[Zeta]]; 
      blim = If[sbound, Quantile[NormalDistribution[0, 1], 1-0.025/maxlag], 
        1.96]; squared = SquaredResiduals /. {opts} /. 
        Options[PlotRacfAR\[Zeta]]; 
      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}], 
        Abs[Sqrt[Transpose[VarianceRacfAR\[Zeta][\[Zeta], lags, maxlag]/n, 
           {1, 1}]]]]; acfrange = {-1, 1}*
        (Ceiling[10*1.06*Max[Abs[r], blim/Sqrt[n]]]/10.); 
      frameticks = Chop[{Range[1, Length[r], Max[1, Round[Length[r]/10.]]], 
         Range[acfrange[[1]], acfrange[[2]], 0.1]}]; 
      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]]];
 
VarianceRacfAR\[Zeta][\[Zeta]_, lags_, m_] := 
  Module[{z, x, y, p = Last[lags], q = Length[lags], X = {}}, 
     z = Array[0 & , p]; 
     Do[z[[lags[[i]]]] = \[Zeta][[i]], {i, q}];
     phi = PacfToAR[z]; 
     x = -Prepend[ARInvert[phi, m - 1], 1]; 
     Do[AppendTo[X, Join[Table[0, {i}], Drop[x, -i]]], {i, 0, p - 1}];
     X = Transpose[X]; 
     Jq = Jacobian\[Phi]\[Zeta][z]; 
     Jq = TakeMatrixColumns[Jq,lags]; 
     XJ = X . Jq; 
     FullSimplify[IdentityMatrix[m]-XJ . Inverse[InformationMatrixZeta[\[Zeta], lags]] . Transpose[XJ]]
];

Jacobian\[Phi]\[Zeta][\[Zeta]_] := 
    Dot @@ Table[\[GothicCapitalJ][\[Zeta], k], {k, Length[\[Zeta]] - 1}] /; 
     Length[\[Zeta]] > 1;
 
Jacobian\[Phi]\[Zeta][\[Zeta]_] := {{1}} /; Length[\[Zeta]] == 1;

(* old version
Jacobian\[Phi]\[Zeta][\[Zeta]_] := 
    Dot @@ Table[\[GothicCapitalJ][\[Zeta], k], {k, Length[\[Zeta]] - 1}];
*)

\[GothicCapitalJ][\[Zeta]_, k_] := 
    Module[{A, J, p = Length[\[Zeta]], \[Phi]pm1j}, 
      J = IdentityMatrix[p - k]; A = ZeroMatrix[p - k, k]; 
       \[Phi]pm1j = Append[FullSimplify[Reverse[
           -PacfToAR[Take[\[Zeta], p - k]]]], 1]; 
       Do[A[[j,1]] = \[Phi]pm1j[[j]], {j, p - k}]; 
       Do[J[[j,p - k + 1 - j]] += -\[Zeta][[p - k + 1]], {j, p - k}]; 
       AppendColumns[AppendRows[J, A], AppendRows[ZeroMatrix[k, p - k], 
         IdentityMatrix[k]]]] /; Length[\[Zeta]] > k && k > 0;
 
\[GothicCapitalJ][\[Zeta]_, k_] := {{1}} /; k == 0;

(* old version 
InfMat\[Zeta][\[Zeta]_] := Module[{J, Ip}, 
     J = Jacobian\[Phi]\[Zeta][\[Zeta]]; 
      Ip = InformationMatrixAR[PacfToAR[\[Zeta]]]; 
      FullSimplify[Transpose[J] . Ip . J]];
 
InfMat\[Zeta][\[Zeta]_, lags_] := Module[{J, Ip, Iz, q = Length[lags], 
      p = Last[lags], z, phi}, z = Array[0 & , p]; 
      FullSimplify[If[p == 1, InformationMatrixAR[\[Zeta]], 
        Do[z[[lags[[i]]]] = \[Zeta][[i]], {i, q}]; 
         J = Jacobian\[Phi]\[Zeta][z]; phi = PacfToAR[z]; 
         Ip = InformationMatrixAR[phi]; Ip = Transpose[J] . Ip . J; 
         Iz = Array[0 & , {q, q}]; Do[Iz[[i,j]] = Ip[[lags[[i]],lags[[j]]]], 
          {i, q}, {j, q}]; Iz]]];
*)

VarianceRacfAR\[Phi][\[Phi]_, lgs_, m_] := Module[{x, y, p, q, lags, X = {}}, 
     lags = Sort[lgs]; p = Last[lags]; q = Length[lags]; 
      phi = Array[0 & , p]; Do[phi[[lags[[i]]]] = \[Phi][[i]], {i, q}]; 
      x = -Prepend[ARInvert[phi, m - 1], 1]; 
      Do[AppendTo[X, Join[Table[0, {i}], Drop[x, -i]]], {i, 0, p - 1}]; 
      X = Transpose[X]; X = TakeMatrixColumns[X, lags]; 
      FullSimplify[IdentityMatrix[m] - X . Inverse[InfMat\[Phi][\[Phi], 
           lags]] . Transpose[X]]];

InfMat\[Phi][\[Phi]_, lg_] := Module[{J, Ip, Iz, p, q, lags}, 
     lags = Sort[lg]; q = Length[lags]; p = Last[lags]; phi = Array[0 & , p]; 
      Do[phi[[lags[[i]]]] = \[Phi][[i]], {i, q}]; 
      Ip = FullSimplify[InformationMatrixAR[phi]]; Iz = Array[0 & , {q, q}]; 
      Do[Iz[[i,j]] = Ip[[lags[[i]],lags[[j]]]], {i, q}, {j, q}]; Iz];

PlotRacfAR\[Phi][res_, \[Phi]_, lags_, (opts___)?OptionQ] := 
    Module[{maxlag, acfrange, frameticks, r, blim, sbound, squared, n}, 
     n = Length[res]; maxlag = MaxLag /. {opts} /. Options[PlotRacfAR\[Phi]]; 
      maxlag = If[maxlag <= Max[lags], 2*Max[lags], maxlag]; 
      sbound = SimultaneousBound /. {opts} /. Options[PlotRacfAR\[Phi]]; 
      blim = If[sbound, Quantile[NormalDistribution[0, 1], 1-0.025/maxlag], 
        1.96]; squared = SquaredResiduals /. {opts} /. 
        Options[PlotRacfAR\[Phi]]; 
      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}], 
        Abs[Sqrt[Transpose[VarianceRacfAR\[Phi][\[Phi], lags, maxlag]/n, 
           {1, 1}]]]]; acfrange = {-1, 1}*
        (Ceiling[10*1.06*Max[Abs[r], blim/Sqrt[n]]]/10.);
      frameticks = Chop[{Range[1, Length[r], Max[1, Round[Length[r]/10.]]], 
         Range[acfrange[[1]], acfrange[[2]], 0.1]}]; 
      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]]
];

(* THEORETICAL AUTOCORRELATION *)

Toeplitz[m_] := Return[Table[m[[1 + Abs[i - j]]], {i, Length[m]}, 
      {j, Length[m]}]];
 
TacvfAR[(phi_List)?NumberVectorQ, MaxLag_Integer] := 
 Module[{p, b, c, rp1}, 
  If[ !And @@ InvertibleQ /@ {phi, {}}, Message[TacvfAR::singular]; Return[{}]]; 
  p = Length[phi]; 
  rp1 = 1 + p; 
  If[p == 0, 
     Prepend[Table[0., {MaxLag}], 1], 
     c = getc[p, 1, phi, {0.}]; 
     b = getb[{0.}, c, 2, rp1]; 
     g = getg[phi, b, p, rp1]; 
     getr[phi, g, p, rp1, MaxLag]
  ]
];

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];
 
InformationMatrixAR[(phi_)?VectorQ] := Module[{vvmatrix, uvmatrix, uumatrix, 
      i, id}, Toeplitz[Drop[Drop[tccf[phi, phi], Length[phi] - 1], -1]]];
 
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[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]]];

AcfPlot[c_, (opts___)?OptionQ] := Module[{acfrange, acfticks, 
      c2 = Prepend[c, 1]}, acfrange = If[MemberQ[(#1 > 0 & ) /@ c2, False], 
        {-1, 1}, {0, 1}]; frameticks = 
       Chop[{Range[0, Length[c2], Max[1, Round[Length[c2]/10.]]], 
         Range[acfrange[[1]], acfrange[[2]], 0.5]}]; 
      Show[Graphics[{RGBColor[0, 0, 1], Thickness[0.01], 
         Line /@ ({{#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, 0.5]}]]];

InformationMatrixZeta[\[Zeta]_] := Module[{J, Ip}, 
     J = If[Length[\[Zeta]]==1, {{1}}, JacobianPhiZeta[\[Zeta]]]; 
     Ip = FullSimplify[InformationMatrixAR[PacfToAR[\[Zeta]]]]; 
     FullSimplify[Transpose[J] . Ip . J]];
 
InformationMatrixZeta[\[Zeta]_, lags_] := 
 Module[{J, Ip, Iz, q = Length[lags], p = Last[lags], z}, 
  z = Array[0 & , p]; 
  Do[z[[lags[[i]]]] = \[Zeta][[i]], {i, q}]; 
      J = If[Length[\[Zeta]]==1, {{1}}, JacobianPhiZeta[\[Zeta]]];
      J = JacobianPhiZeta[\[Zeta]];
      Print["J=",J];
      Ip = FullSimplify[InformationMatrixAR[
         PacfToAR[z]]]; Ip = FullSimplify[Transpose[J] . Ip . J]; 
      Iz = Array[0 & , {q, q}]; Do[Iz[[i,j]] = Ip[[lags[[i]],lags[[j]]]], 
       {i, q}, {j, q}]; Iz];

InformationMatrixZeta[\[Zeta]_, lags_] := 
 Module[{J, Ip, Iz, q = Length[lags], p = Last[lags], z}, 
  z = Array[0 & , p]; 
  Do[z[[lags[[i]]]] = \[Zeta][[i]], {i, q}]; 
 J = JacobianPhiZeta[z]; 
 J = If[Length[z]==1, {{1}}, JacobianPhiZeta[z]];
 Ip = FullSimplify[InformationMatrixAR[PacfToAR[z]]]; 
 Ip = FullSimplify[Transpose[J] . Ip . J]; 
 Iz = Array[0 & , {q, q}]; 
 Do[Iz[[i,j]] = Ip[[lags[[i]],lags[[j]]]],{i, q}, {j, q}]; 
Iz];
 
 
JacobianPhiZeta[\[Zeta]_] := Dot @@ Table[JacobianK[\[Zeta], k], 
      {k, Length[\[Zeta]] - 1}];
 
JacobianK[\[Zeta]_, k_] := Module[{A, J, p = Length[\[Zeta]], \[Phi]pm1j}, 
      J = IdentityMatrix[p - k]; A = ZeroMatrix[p - k, k]; 
       \[Phi]pm1j = Append[FullSimplify[Reverse[
           -PacfToAR[Take[\[Zeta], p - k]]]], 1]; 
       Do[A[[j,1]] = \[Phi]pm1j[[j]], {j, p - k}]; 
       Do[J[[j,p - k + 1 - j]] += -\[Zeta][[p - k + 1]], {j, p - k}]; 
       AppendColumns[AppendRows[J, A], AppendRows[ZeroMatrix[k, p - k], 
         IdentityMatrix[k]]]] /; Length[\[Zeta]] > k;
   
(* SPECTRAL DENSITY FUNCTION *)

SdfARPlot[(phi_)?NumberVectorQ, (opts___)?OptionQ] := 
    Module[{pr, pl, sdf, framelab, xlog, ylog, d, f, n}, 
     xlog = LoggedF /. {opts} /. Options[SdfARPlot]; 
      ylog = LoggedSdf /. {opts} /. Options[SdfARPlot]; 
      v = InnovationVariance /. {opts} /. Options[SdfARPlot]; 
      pr = PlotRange /. {opts} /. Options[SdfARPlot]; 
      pl = PlotLabel /. {opts} /. Options[SdfARPlot]; 
      RFQ = RadialFrequency /. {opts} /. Options[SdfARPlot]; 
      {f, sdf} = Transpose[SdfAR[phi, InnovationVariance -> v]];
      f = If[RFQ, 2 Pi f, f]; 
      f = If[xlog, Log[2, f], f]; 
      sdf = If[ylog, Log[2, sdf], sdf]; 
      framelab = Which[ 
                 !(xlog || ylog), {"\[Lambda]", "f(\[Lambda])"}, 
                 !xlog && ylog, {"\[Lambda]", "\!\(log\_2\)f(\[Lambda])"}, 
		 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 -> True]];

SdfAR[(phi_)?NumberVectorQ, (opts___)?OptionQ] := 
    Module[{pFFT = 9, v, \[ScriptP]}, 
     v = InnovationVariance /. {opts} /. Options[SdfAR]; 
      \[ScriptP] = v/(2^(pFFT + 1)*Re[#1*Conjugate[#1]] & )[
         Take[Fourier[Join[{1}, -phi, Table[0, {2^(pFFT + 1) - 1 - Length[
                phi]}]]], 2^pFFT]]; Transpose[{Range[0, 511]/1024., 
        \[ScriptP]}]];
 
SdfAR[(phi_)?VectorQ, \[Lambda]_, (opts___)?OptionQ] := 
    Module[{v, norq}, v = InnovationVariance /. {opts} /. Options[SdfAR]; 
      v/If[Length[phi] > 0, Abs[1 - Sum[phi[[k]]*Exp[-2*Pi*I*\[Lambda]*k], 
            {k, Length[phi]}]], 1]^2];

(* SAMPLE AUTOCORRELATION *)

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, 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.]]], 
        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]]}, {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]]];
 
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]
    ];
 
(* PARTIAL AUTOCORRELATION *)

SPacfPlot[(z_)?NumberVectorQ, (opts___)?OptionQ] := 
  Module[{slim,sbound,method, ifmean0, maxlag, cilim, frameticks, r}, 
   method = Method /. {opts} /. Options[SPacfPlot]; 
    ifmean0 = MeanZero /. {opts} /. Options[SPacfPlot]; 
    maxlag = MaxLag /. {opts} /. Options[SPacfPlot]; 
    sbound = SimultaneousBound /. {opts} /. Options[SPacfPlot]; 
    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", "SPACF"}, 
      PlotRange -> {-1, 1}, FrameTicks -> frameticks, 
      GridLines -> {{}, {0.}}]]
      ];



Jacobian\[Phi]\[Zeta][\[Zeta]_] := 
    Dot @@ Table[\[GothicCapitalJ][\[Zeta], k], {k, Length[\[Zeta]] - 1}];
 
\[GothicCapitalJ][\[Zeta]_, k_] := 
    Module[{A, J, p = Length[\[Zeta]], \[Phi]pm1j}, 
      J = IdentityMatrix[p - k]; A = ZeroMatrix[p - k, k]; 
       \[Phi]pm1j = Append[FullSimplify[Reverse[
           -PacfToAR[Take[\[Zeta], p - k]]]], 1]; 
       Do[A[[j,1]] = \[Phi]pm1j[[j]], {j, p - k}]; 
       Do[J[[j,p - k + 1 - j]] += -\[Zeta][[p - k + 1]], {j, p - k}]; 
       AppendColumns[AppendRows[J, A], AppendRows[ZeroMatrix[k, p - k], 
         IdentityMatrix[k]]]] /; Length[\[Zeta]] > k;
      

PacfPlot[(z_)?NumberVectorQ, (opts___)?OptionQ] := 
    Module[{ifmean0, maxlag, cilim, frameticks, nr, r, rsd, rlo, gop, g, 
      LargeP, OtherP, Lr, lgest, LP, OP, LargePSorted, ycor, MLEQ, YWQ}, 
     ifmean0 = MeanZero /. {opts} /. Options[PacfPlot]; 
      maxlag = MaxLag /. {opts} /. Options[PacfPlot]; 
      npanels = NumberOfPanels /. {opts} /. Options[PacfPlot]; 
      lgest = Largest /. {opts} /. Options[PacfPlot]; 
      indmthd = Method /. {opts} /. Options[PacfPlot]; 
      MLEQ = Which[indmthd === Automatic && maxlag <= 30, True, indmthd === "MLE", 
        True, True, False]; YWQ = If[indmthd === "YW", True, False]; 
      nr = Floor[maxlag/npanels]; ycor = If[ifmean0, z, z - Mean[z]]; 
      r = Which[MLEQ, Last[GetARFit[Range[maxlag], ycor]], YWQ, 
        PacfDL[Sacvf[ycor, MaxLag -> maxlag]], True, PacfBurg[ycor, maxlag]];
      rsd = Sqrt[Transpose[Inverse[Length[z]*InformationMatrixZeta[r]], {1, 1}]]; 
      rlo = r - 1.96*rsd; rhi = r + 1.96*rsd; 
      LargeP = Take[Sort[Transpose[{Range[Length[r]], Abs[r]}], 
         Last[#1] > Last[#2] & ], lgest]; LargeP = If[LargeP == {}, {}, 
        First[Transpose[LargeP]]]; LargePSorted = Sort[LargeP]; 
      OtherP = Complement[Range[maxlag], LargePSorted]; 
      gop = {Frame -> True, Axes -> False, AxesOrigin -> {0, 0}, 
        FrameLabel -> {"Lag", "Acf"}, PlotRange -> {-1, 1}, 
        FrameTicks -> {Range[0, Length[r], Max[1, Round[nr/10.]]], 
          Range[-1, 1, 0.5]}, GridLines -> {{}, Range[-1, 1, 0.5]}}; 
      g = If[npanels <= 1, Graphics[{{RGBColor[0, 0, 1], Thickness[0.01], 
           Line /@ ({{#1, rlo[[#1]]}, {#1, rhi[[#1]]}} & ) /@ OtherP}, 
          {RGBColor[1, 0, 0], Thickness[0.01], Line /@ 
            ({{#1, rlo[[#1]]}, {#1, rhi[[#1]]}} & ) /@ LargePSorted}}, gop], 
        Lr = Partition[Range[maxlag], nr]; 
         OP = Table[Complement[Lr[[i]], LargePSorted], {i, Length[Lr]}]; 
         LP = Table[Intersection[Lr[[i]], LargePSorted], {i, Length[Lr]}]; 
         GraphicsArray[Partition[Table[Graphics[{{RGBColor[0, 0, 1], 
               Thickness[0.01], Line /@ ({{#1, rlo[[#1]]}, {#1, rhi[[
                     #1]]}} & ) /@ OP[[i]]}, {RGBColor[1, 0, 0], Thickness[
                0.01], Line /@ ({{#1, rlo[[#1]]}, {#1, rhi[[#1]]}} & ) /@ 
                 LP[[i]]}}, gop], {i, npanels}], 1]]]; Show[g]; 
      If[lgest == 0, g, {g, LargeP}]];
 

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;

ICFromPacf[n_Integer, (pacf_)?NumberVectorQ, (opts___)?OptionQ] := 
    Module[{tableQ, indtable, k, ic, whichic, neg2LL, plotQ}, 
      whichic = Criterion /. {opts} /. Options[ICFromPacf]; 
       indtable = Table /. {opts} /. Options[ICFromPacf]; 
       tableQ = If[indtable === False, False, True]; 
       k = If[IntegerQ[indtable] && indtable > 0, indtable, 5]; 
       plotQ = Plot /. {opts} /. Options[ICFromPacf]; 
       neg2LL = DetAR[pacf] + 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]]]; 
       If[plotQ, ListPlot[Transpose[{Range[0, Length[pacf]], ic}], 
         Axes -> False, Frame -> True, PlotRange -> All, 
         PlotStyle -> {AbsolutePointSize[6]}]]; 
       If[tableQ, Take[Sort[Transpose[{Range[0, Length[ic] - 1], ic}], 
          Last[#1] < Last[#2] & ], k], IndexOfMinimum[ic] - 1]] /; n > 0;
 
PacfToAR[(pi_)?VectorQ /; Length[pi] > 0] := 
    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];
 
PacfToAR[pi_] := {};

ARToPacf[(phi_)?VectorQ] := Module[{pi, a, phikp1, phik, L}, 
     swQ = False; phik = phi; L = Length[phi]; pi = Table[0, {L}]; 
      Do[pi[[L + 1 - k]] = a = Last[phik]; phikp1 = Drop[phik, -1]; 
        swQ = Xor[swQ, Abs[a] === 1]; phik = FullSimplify[
          If[swQ, Array[0 & , Length[phikp1]], (phikp1 + a*Reverse[phikp1])/
            (1 - a^2)]]; , {k, 1, L}]; pi];

PacfDL[(c_)?NumericVectorQ, (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]];

AICSubAR[z_, (opts___)?OptionQ] := Module[{n = Length[z], indmQ, maxlag, K, 
      r, LargeP, ac, bc, p, LPi, OtherP, ri, neg2LL, grph1, grph2, tb, 
      minAIC, minBIC, indmthd, MLEQ, YWQ}, 
     indmQ = MeanZero /. {opts} /. Options[AICSubAR]; 
      maxlag = MaxLag /. {opts} /. Options[AICSubAR]; 
      maxlag = If[maxlag === Automatic, Min[Ceiling[n/4], 200], maxlag]; 
      K = Largest /. {opts} /. Options[AICSubAR]; 
      indmthd = Method /. {opts} /. Options[AICSubAR]; 
      MLEQ = Which[indmthd === Automatic && maxlag <= 5, True, 
        indmthd === "MLE", True, True, False]; 
      YWQ = If[indmthd === "YW", True, False]; 
      ycor = If[indmQ, z, z - Mean[z]]; 
      r = Which[MLEQ, Last[GetARFit[Range[maxlag], ycor]], YWQ, 
        PacfDL[Sacvf[ycor, MaxLag -> maxlag]], True, PacfBurg[ycor, maxlag]]; 
      LargeP = Take[Sort[Transpose[{Range[Length[r]], Abs[r]}], 
         Last[#1] > Last[#2] & ], K]; LargeP = If[LargeP == {}, {}, 
        First[Transpose[LargeP]]]; ac = bc = {0}; p = {{0}}; 
      Do[LPi = Take[LargeP, i]; p = Append[p, LPi]; 
        OtherP = Complement[Range[maxlag], LPi]; ri = r; ri[[OtherP]] = 0; 
        ac = {ac, AICFromPacf[ri, n, i]}; bc = {bc, BICFromPacf[ri, n, i]}; , 
       {i, K}]; ac = Flatten[ac]; bc = Flatten[bc]; ac = Rest[ac]; 
      bc = Rest[bc]; p = Rest[p]; grph1 = ListPlot[Transpose[{Range[K], ac}], 
        Axes -> False, Frame -> True, PlotRange -> All, 
        PlotStyle -> {AbsolutePointSize[6]}, FrameLabel -> {"k", "AIC"}]; 
      grph2 = ListPlot[Transpose[{Range[K], bc}], Axes -> False, 
        Frame -> True, PlotRange -> All, PlotStyle -> {AbsolutePointSize[6]}, 
        FrameLabel -> {"k", "BIC"}]; 
      tb = TableForm[MapThread[Append, {Transpose[{Range[K], ac, bc}], p}], 
        TableHeadings -> {None, {"k", "AIC", "BIC", "Model"}}, 
        TableDepth -> 2]; minAIC = p[[IndexOfMinimum[ac]]]; 
      minBIC = p[[IndexOfMinimum[bc]]]; {grph1, grph2, SubsetTable -> tb, 
       MinAICModel -> minAIC, MinBICModel -> minBIC}];


(* FORECAST *)

ForecastOneStepTestData[ans_, yTest_] := 
    Module[{\[Phi]s, mu, z, zTest, p, zX, X, zP, err, MSEerr, MSEtot, tz, tb, 
      eps, rz, g}, \[Phi]s = ARCoefficients /. ans; mu = MeanValue /. ans; 
      z = (OriginalSeries /. ans) - mu; zTest = yTest - mu; 
      nT = Length[zTest]; p = Length[\[Phi]s]; 
      zX = Flatten[{Take[Reverse[z], p], Drop[zTest, -1]}]; 
      X = Reverse /@ Table[Take[Drop[zX, i], p], {i, 0, nT - 1}]; 
      zP = X . \[Phi]s; err = zTest - zP; MSEerr = Total[err^2]/Length[err]; 
      MSEtot = Total[zTest^2]/Length[zTest]; zTest = zTest + mu; 
      zP = zP + mu; tz = Length[z] + Range[nT] - 1; 
      tb = TableForm[Transpose[{zTest, zP, err}], TableHeadings -> 
         {tz, {"Obs.", "Fit", "Error"}}]; 
      eps = 0.04*(Max[zTest, zP] - Min[zTest, zP]); 
      rz = {Min[zTest, zP] - eps, Max[zTest, zP] + eps}; 
      g = ListPlot[Transpose[{tz, zTest}], PlotJoined -> True, 
        PlotRange -> rz, Axes -> False, Frame -> True, 
        FrameLabel -> {"observation number", 
          "data (line), forecast (points)"}, Epilog -> {AbsolutePointSize[5], 
          RGBColor[0, 0, 1], Point /@ Transpose[{tz, zP}]}]; 
      {g, ForecastTable -> tb, RMSE -> Sqrt[MSEerr], 
       RSQ -> Round[((MSEtot - MSEerr)/MSEtot)*100], XMatrix -> X}];
  
(* SIMULATION *)

SimulateGaussianAR[(phi_)?InvertibleQ, n_Integer/;n>0]:=
    Module[
      {p=Length[phi],a,zp},
      a=RandomArray[NormalDistribution[0,1],n-p];
      If[p==0,Return[a]];
      zp=Random[
          MultinormalDistribution[Table[0,{p}],Toeplitz[TacvfAR[phi,p-1]]]];
      FastGetSim[zp,a,phi]
      ];
 
SimulateGaussianAR[phi_?NumericQ, n_Integer /; n > 0] := 
    SimulateGaussianAR[{phi}, n];
 
SimulateGaussianAR[ans_] := (MeanValue /. ans) + (ResidualSd /. ans)*
      SimulateGaussianAR[ARCoefficients /. ans, Length[OriginalSeries /. ans]];

BootAR[ans_] := Module[{sd, zm, phi, n}, sd = ResidualSd /. ans; 
      zm = MeanValue /. ans; phi = ARCoefficients /. ans; 
      n = Length[OriginalSeries /. ans]; zm + SimulateGaussianAR[phi, n]*sd];
 
 
(* UTILITIES *)
ToeplitzQ[x_] := SquareMatrixQ[x] && x == Toeplitz[x[[1]]];

SquareMatrixQ[x_] := MatrixQ[x] && Apply[Equal, Dimensions[x]];

AICFromPacf[pacf_, n_, k_] := n*Log[Times @@ (1 - pacf^2)] + 2*k;
 
BICFromPacf[pacf_, n_, k_] := n*Log[Times @@ (1 - pacf^2)] + Log[n]*k;

NumberVectorQ[x_] := VectorQ[x] && And @@ NumberQ /@ x;

NumericVectorQ[x_] := VectorQ[x] && And @@ NumericQ /@ x;

TakeMatrixColumns[j_, lags_] := Module[{i = Sort[lags], q = Length[lags], jt, 
      jt2 = {}, p = First[Dimensions[j]]}, jt = Transpose[j]; 
      Do[jt2 = {jt2, Take[jt, {i[[k]]}]}, {k, q}]; 
      Transpose[Partition[Flatten[jt2], p]]];

IndexOfMinimum[(x_)?NumberVectorQ] := 
    First[First[Select[Transpose[{Range[Length[x]], x}], 
       #1[[2]] == Min[x] & ]]];

InvertibleQ[phi_] := With[{Fuzz = 10^(-14)}, 
     NumericVectorQ[phi]&&(And @@ (Abs[#1] + Fuzz < 1. & ) /@ ARToPacf[phi])];
   
ARInvert[(phi_)?VectorQ, MaxLag_ /; MaxLag > 0] := 
    Module[{B}, Take[Drop[CoefficientList[
        Normal[Series[1/ARPolynomial[phi, B], {B, 0, 2*MaxLag}]], B], 1], 
      MaxLag]];
 
ARPolynomial[phi_, B_] := Prepend[-phi, 1] . Table[B^j, {j, 0, Length[phi]}];

ARRoots[phi_] := Which[Length[phi] > 1, 
     Table[NRoots[1 - Sum[phi[[i]]*B^i, {i, 1, Length[phi]}] == 0, B][[i,2]], 
      {i, 1, Length[phi]}], Length[phi] == 1, {1/First[phi]}, True, {}];

ReadData[(datafile_)?StringQ, (opts___)?OptionQ] := 
    Module[{continuetitle = True, z = {}, i, c, title, zt, NoTitle = True}, 
     If[ !File === FileType[datafile], Message[ReadData::NotFound, datafile]; 
        Return[{}]]; sQ = Silent /. {opts} /. Options[ReadData]; 
      OpenRead[datafile]; title = ""; While[continuetitle, 
       i = StreamPosition[datafile]; c = ToString[Read[datafile, Character]]; 
        c = Which[LetterQ[c], c, c == "#", "", True, 0]; 
        If[c != 0, NoTitle = False]; If[c == 0, continuetitle = False, 
         title = StringJoin[StringJoin[title, "\n"], Read[datafile, String, 
            RecordSeparators -> "\n"]]]]; SetStreamPosition[datafile, i]; 
      While[ToString[zt = Read[datafile, Real]] != "EndOfFile", 
       AppendTo[z, zt]]; Close[datafile]; 
      If[ !sQ, If[ !NoTitle, Print[StringDrop[title, 1]]]; 
        Print["Number of Values Input = ", Length[z]]]; z];
 
(* Compiled Functions *)

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\[Equal]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]]]]];

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]]]];

Qd=Compile[{{phi,_Real,1},{ChampD,_Real,2}},phi.ChampD.phi];

Dta = Compile[{{z, _Real, 1}},Times@@((1 - z^2)^-Range[Length[z]])];

DtaFUN[x_] := Dta[x]; (* do not change this! *)

FastGetSim=
    Compile[{{zp, _Real, 1}, {a, _Real,1},{phi, _Real, 1}},
      Join[Reverse[zp],
        Rest[FoldList[Prepend[Drop[#1,-1],phi.#1+#2]&,zp,a]]//Transpose//
          First]] ;

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

End[];
EndPackage[];
