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

(* Author:  *)
(* Copyright (c), A.I. McLeod, 2005 *)
(* ARMA is a package for fitting arma models *)
(* Summary:
*)


Needs["Statistics`DescriptiveStatistics`"];

(********************************************************************) 
BeginPackage["FitARMA`", "FitAR`"]
(********************************************************************) 

(* global variables *)
ARMAModel; Loglikelihood; ParameterTable; ARMACoefficients; MeanValue; 
ResidualSd; PortmanteauDiagnostic; AdjustableParameters; Residuals; OriginalSeries;
ARMAModelOrder; 
NumIterations;Mu;ErrorTolerance;

(* FitAR.m must be placed in the Mathematica Applications Directory for this to work *)
(* Needs["FitAR`"] *)
 
(* Load needed Standard Packages *)
Needs["Statistics`DescriptiveStatistics`"];
Needs["Statistics`ContinuousDistributions`"];
Needs["Statistics`DataManipulation`"];
Needs["Statistics`MultinormalDistribution`"]; 
Needs["LinearAlgebra`MatrixManipulation`"];

(* FitARMA *)

Options[FitARMA] := {MaxLag -> Automatic, MeanZero -> False, 
     ARApproximation -> 30};
 
FitARMA /: FitARMA::Invalidpdq = 
     "Input pdq must be either a vector of length 2, 3 or 7 or else a \
number.";
 
FitARMA /: FitARMA::usage = "ans=FitARMA[z, pdq] fits an ARIMA(p,d,q) to \
series z, where pdq={p,d,q}. To fit a \
SARIMA(p,d,q)(\!\(p\_s\),\!\(d\_s\),\!\(q\_s\)\!\(\()\_s\)\) use \
pdq={p,d,q,\!\(p\_s\),\!\(d\_s\),\!\(q\_s\),s}. Also for convenience you may \
use pdq=p for AR(p) and pdq={p,q} for ARMA. The result of FitARIMA, ans, is a \
lengthy rule containing the output. It is best to use one or more of the \
following auxiliary functions to extract information from ans. Summary[ans] \
terse summary of the fitted model, Coef[ans] gives a table of the parameter \
coefficients, DiagnosticPlotARIMA[ans] produces one or more diagnostic plots. \
See the help function for each of these functions for more information.";
 
Options[GetARMAFit] = {InitialEstimates -> 0};

GetARMAFit /: GetARMAFit::switch = "Switch to RandomSearch";
 
GetARMAFit /: GetARMAFit::usage = "GetARMAFit[p,q,pApprox,y] fits an \
ARMA(p,q) to mean-corrected data y using high order AR(pApprox) \
approximation. Note: it returns the transformed parameters so to get the \
usual ARMA parameters you need to use PacfToAR to transform back. Takes \
option InitialEstimates which has default value 0.  Otherwise if set to a \
vector of p+q these are used as initial estimates. The initial estimates are \
assumed to be transformed to be in the range (-1,1). In particular, \
InitialEstimates->HRARMA[p,q,y] will use Hannan-Riassen estiamtes.";


Options[GetARMAMeanFit] = {InitialEstimates -> 0};
 
GetARMAMeanFit /: GetARMAMeanFit::usage = "GetARMAMeanFit[p,q,pApprox,z] fits \
ARMA(p,q) including mean parameter, mu.  Uses maximim likelihood estimation \
for all parameters.  Returns a list of the form {loglikelihood, {phi,theta}, \
Mu->mean, NumIterations->n, ErrorTolerance->err.  Takes option \
InitialEstimates which has default value 0. Otherwise it is the initial \
estimates for the phi and theta parameters in the transformed domain -1 to \
+1.  As for example the initial estimates obtained from HRARMA.";

Options[HRARMA] = {AROrder -> 20};
 
HRARMA /: HRARMA::usage = "HRARMA[p,q,x] asymptotically efficient estimates of \
the phi and theta parameters in an ARMA(p,q) model by fitting a high-order AR \
model using the Burg algorithm. The method of Hannan and Rissanen (1982) is \
used. Takes option AROrder which as default value of 20.";
 
getAR\[Mu] /: getAR\[Mu]::usage = "getAR\[Mu][z,phi] computes the exact mle \
for the  mean given parameters phi in an AR model.";
 
(* THEORETICAL AUTOCORRELATION and INFORMATION MATRIX *)

Options[TacvfARMA] = {FractionalDifference -> 0, LagTruncation -> 512};
 
TacvfARMA /: TacvfARMA::singular = "Non-invertible AR or MA model or very \
close to it: the A Matrix is singular or ill-conditioned.";
 
TacvfARMA /: TacvfARMA::unavail = 
     "Symbolic fractional autocovariances are not available";
 
TacvfARMA /: 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.";

(* INFORMATION MATRIX *)

 
(* SPECTRAL DENSITY FUNCTION *)

SdfARut::usage="Utility fn. SdfAR[phi]";

SdfMAut::usage="Utility fn. SdfMA[theta]";

SdfFDWNut::usage="Utility fn. SdfFDWN[d,n]";

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


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

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

(* SIMULATION *)

Options[SimulateARIMA] = {Differencing -> 0, InnovationVariance -> 1};
 
SimulateARIMA /: 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 \
option Differencing with default value 0. 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.";
 
Options[SimulateGaussianARMA] = {InnovationVariance -> 1};
 
SimulateGaussianARMA /: SimulateGaussianARMA::usage = "SimulateGaussianARMA[p\
hi, 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. The innovation variance is assumed 1 but more generally you \
can use SimulateGaussianARMA[phi,theta,n,InnovationVariance->v]";

SimulateNonGaussianARMA /: SimulateNonGaussianARMA::usage = "SimulateNonGauss\
ianARMA[phi,theta,Q,a] simulates a non-Gaussian ARMA with parameters phi and \
theta which must be vectors. Q is the order of the MA approximation. a is the \
vector of innovations of length n+Q, where n is the length of time series \
output.";

SimMA /: SimMA::usage = "SimMA[psi, a] given the vector of length Q+1 of \
impulse coefficients {1, psi_1, ..., psi_Q} and a vector a of length n+Q of \
innovations, the MA(Q) model is simulated.  The computation uses the \
ListConvolve function.";

(* INFORMATION MATRIX *)

InformationMatrixARMA /: InformationMatrixARMA::usage = "InformationMatrixARM\
A[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.";
 

(* DIAGNOSTIC CHECKS AND SUMMARIES *)

CoefARMA /: CoefARMA::usage = "CoefARMA[ans], where ans is the output from \
FitARMA, produces a table showing the fitted model parameters.";
 
SummaryARMA /: SummaryARMA::usage = "SummaryARMA[ans], where ans is the output from \
FitAR, produces a summary of the fitted model.";

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


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

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

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

(*=========================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.";
 
(* UTILITIES *)
BoundaryTestQ /: BoundaryTestQ::usage = "BoundaryTestQ[phi] tests if phi \
parameters are on the boundary. Returns True if there is at least one \
parameter on the boundary. Symbolic values are not allowed. Takes argument \
eps with default value 0.";
 
Options[InvertedCoefficientsARMA] = {Difference -> 0, SeasonalPeriod -> 1, 
     SeasonalDifference -> 0};
 
InvertedCoefficientsARMA /: InvertedCoefficientsARMA::usage = "InvertedCoeffi\
cientsARMA[phi, theta, MaxLag] returns the vector  \
{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.";

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

ImpulseCoefficientsARMA::usage=
"ImpulseCoefficientsARMA[phi,theta,m] yields {1, psi_1, psi_2, ...} where z_t = a_t + psi_1 a_{t-1} + ...";

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

DetAR /: DetAR::usage = 
     "DetAR[phi] computes the normalized covariance determinant of AR.";

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

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

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

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

(* COMPILED FUNCTIONS *)

Qd /: Qd::usage = "Compiled utility function"



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

(*FITARMA AND RELATED FUNCTIONS*)

FitARMA[(y_)?NumberVectorQ, (pdq_)?NumberVectorQ, (opts___)?OptionQ] := 
    Module[{}, ml = MaxLag /. {opts} /. Options[FitARMA]; 
      indmean = MeanZero /. {opts} /. Options[FitARMA]; 
      pApprox = ARApproximation /. {opts} /. Options[FitARMA]; 
      {p, d, q, ps, ds, qs, s} = Which[IntegerQ[psd] && pdq >= 0, 
        {pdq, 0, 0, 0, 0, 0, 0}, Length[pdq] == 1, {First[pdq], 0, 0, 0, 0, 
         0, 0}, Length[pdq] == 2, {pdq[[1]], 0, pdq[[2]], 0, 0, 0, 0}, 
        Length[pdq] == 3, Join[pdq, {0, 0, 0, 0}], True, 
        Message[FitARMA::Invalidpdq]; Return[{}]]; n = Length[y]; 
      ml = If[ml === Automatic, Which[n <= 100, 15, n <= 200, 20, n <= 300, 
         25, n <= 500, 40, True, 50], ml]; mu = If[indmean, 0, Mean[y]]; 
      ycor = y - mu; ans = GetARMAFit[p, q, ps, qs, s, pApprox, ycor]; 
      LL = First[ans]; {ztaPhi, ztaTheta, ztaSeasonalPhi, ztaSeasonalTheta} = 
       If[s > 0, Last[ans], Join[Last[ans], {{}, {}}]]; 
      {phis, thetas, Seasonalphis, Seasonalthetas} = 
       PacfToAR /@ {ztaPhi, ztaTheta, ztaSeasonalPhi, ztaSeasonalTheta}; 
      \[DoubleStruckCapitalD] = GetARMA\[DoubleStruckCapitalD][p, q, ps, qs, 
        s, pApprox, ycor]; resSd = ResidualSDARMA[n, \[DoubleStruckCapitalD], 
        ztaPhi, ztaTheta, ztaSeasonalPhi, ztaSeasonalTheta, s]; 
      Var\[Beta] = InformationMatrixARMA[phis, thetas, Seasonalphis, 
        Seasonalthetas]; se\[Beta] = Sqrt[Transpose[Inverse[Var\[Beta]]/n, 
         {1, 1}]]; sdmn = resSd*Sqrt[1/n]*
        Sqrt[(sumsqPhi[thetas]*sumsqPhi[Seasonalthetas])/
          (sumsqPhi[phis]*sumsqPhi[Seasonalphis])]; 
      se = Append[se\[Beta], sdmn]; 
      beta = Flatten[{{phis, thetas, Seasonalphis, Seasonalthetas}, mu, 
         resSd}]; trows = Flatten[{Table[Subscript["\[Phi]", i], {i, p}], 
         Table[Subscript["\[Theta]", i], {i, q}], 
         Table[Subscript["\[CapitalPhi]", i], {i, ps}], 
         Table[Subscript["\[CapitalTheta]", i], {i, qs}], 
         {"\[Mu]", "\!\(\[Sigma]\_a\)"}}]; 
      tb = TableForm[Transpose[{beta, Append[se, " "], 
          Append[Drop[beta, -1]/se, " "]}], TableHeadings -> 
         {trows, {"Parameter", "SE", "Z"}}]; LL = First[ans]; 
      k = p + q + ps + qs + 1; ai = -2*LL + 2*k; bi = -2*LL + Log[n]*k; 
      res = BackcastResiduals[y - mu, Sequence @@ SARMAToARMA[phis, thetas, 
          Seasonalphis, Seasonalthetas, s]]; 
      pt = PortmanteauTest[res, ml, k - 1, TestType -> "BoxLjung"]; 
      {ARMAModel -> {p, d, q, ps, ds, qs, s}, Loglikelihood -> LL, 
       ParameterTable -> tb, ARMACoefficients -> Drop[beta, -2], 
       MeanValue -> mu, ResidualSd -> resSd, PortmanteauDiagnostic -> pt, 
       AdjustableParameters -> k, Residuals -> res, OriginalSeries -> y}];
 
 
GetARMAFit[(p_)?NumberQ, (q_)?NumberQ, (ps_)?NumberQ, (qs_)?NumberQ, 
     (s_)?NumberQ, (pApprox_)?NumberQ, (y_)?NumberVectorQ] := 
    Module[{\[DoubleStruckCapitalD], zinit, zeps, vars, n = Length[y]}, 
     If[ps == 0 && qs == 0, GetARMAFit[p, q, pApprox, y], 
      \[DoubleStruckCapitalD] = If[q == 0 && qs == 0, ChampernowneD[y, 
          p + ps, MeanZero -> True], ChampernowneD[y, papprox, 
          MeanZero -> True]]; ztaAR = Table[Subscript[\[ScriptZ], i], 
         {i, p}]; ztaMA = Table[Subscript[\[ScriptZ], i + p], {i, q}]; 
       ztasAR = Table[Subscript[\[ScriptZ], i + p + q], {i, ps}]; 
       ztasMA = Table[Subscript[\[ScriptZ], i + p + q + ps], {i, qs}]; 
       zinit = Array[0. & , p + q + ps + qs]; zeps = Min[1 - Abs[zinit]]*0.5; 
       vars = Table[{Subscript[\[ScriptZ], i], -zeps + zinit[[i]], 
          zeps + zinit[[i]], -0.99999, 0.99999}, {i, p + q + ps + qs}]; 
       ans = NMaximize[{FastLoglikelihoodARMA[Length[y], 
           \[DoubleStruckCapitalD], ztaAR, ztaMA, ztasAR, ztasMA], 
          Evaluate[Sequence @@ cons]}, Evaluate[v], Method -> "NelderMead"]; 
       {First[ans], {ztaAR, ztaMA, ztasAR, ztasMA} /. Last[ans]}]];

GetARMAFit[(p_)?NumberQ, (q_)?NumberQ, (pApprox_)?NumberQ, 
     (y_)?NumberVectorQ, (opts___)?OptionQ] := 
    Module[{\[DoubleStruckCapitalD], init, zinit, zeps, vars, n = Length[y], 
      ztaAR, ztaMA, ans}, init = InitialEstimates /. {} /. 
        Options[GetARMAFit]; zinit = If[init === 0, Array[0. & , p + q], 
        init]; \[DoubleStruckCapitalD] = Which[p == 0 && q == 0, {{y . y}}, 
        q == 0, ChampernowneD[y, p, MeanZero -> True], q > 0, 
        ChampernowneD[y, pApprox, MeanZero -> True]]; 
      If[p == 0 && q == 0, {FastLoglikelihoodARMA[n, \[DoubleStruckCapitalD], 
         {}, {}], {}, {}}, ztaAR = Table[Subscript[\[ScriptZ], i], {i, p}]; 
        ztaMA = Table[Subscript[\[ScriptZ], i + p], {i, q}]; 
        zeps = Min[1 - Abs[zinit]]*0.5; 
        vars = Table[{Subscript[\[ScriptZ], i], -zeps + zinit[[i]], 
           zeps + zinit[[i]], -0.99999, 0.99999}, {i, p + q}]; 
        ans = Check[FindMaximum[Hold[FastLoglikelihoodARMA[Length[y], 
             \[DoubleStruckCapitalD], ztaAR, ztaMA]], 
           Evaluate[Sequence @@ vars], Compiled -> False], 
          cons = Table[-1 < Subscript[\[ScriptZ], i] < 1, {i, p + q}]; 
           v = Table[Subscript[\[ScriptZ], i], {i, p + q}]; 
           NMaximize[{FastLoglikelihoodARMA[Length[y], 
              \[DoubleStruckCapitalD], ztaAR, ztaMA], Evaluate[
              Sequence @@ cons]}, Evaluate[v], Method -> "NelderMead"], 
          FindMaximum::regex, FindMaximum::regex1]; 
        ans = If[NumericQ[First[ans]], ans, Message[GetARMAFit::switch]; 
           cons = Table[-1 < Subscript[\[ScriptZ], i] < 1, {i, p + q}]; 
           v = Table[Subscript[\[ScriptZ], i], {i, p + q}]; 
           NMaximize[{FastLoglikelihoodARMA[Length[y], 
              \[DoubleStruckCapitalD], ztaAR, ztaMA], Evaluate[
              Sequence @@ cons]}, Evaluate[v], Method -> "RandomSearch"]]; 
        {First[ans], {ztaAR, ztaMA} /. Last[ans]}]];
 
FastLoglikelihoodARMA[n_, \[DoubleStruckCapitalD]_, (ztaAR_)?NumberVectorQ, 
     (ztaMA_)?NumberVectorQ] := Module[{phiApprox, nl, papprox}, 
     -Which[Length[ztaMA] == 0 && Length[ztaAR] == 0, 
       0.5*n*Log[\[DoubleStruckCapitalD][[1,1]]/n], Length[ztaMA] == 0, 
       0.5*Log[Dta[ztaAR]] + 0.5*n*Log[Qd[Prepend[PacfToAR[ztaAR], -1.], 
            \[DoubleStruckCapitalD]]/n], Length[ztaAR] == 0, 
       papprox = First[Dimensions[\[DoubleStruckCapitalD]]] - 1; 
        phiApprox = InvertedCoefficientsARMA[{}, PacfToAR[ztaMA], papprox]; 
        0.5*Log[Dta[ztaMA]] + 0.5*n*Log[Qd[Prepend[phiApprox, -1.], 
             \[DoubleStruckCapitalD]]/n], True, 
       papprox = First[Dimensions[\[DoubleStruckCapitalD]]] - 1; 
        phiApprox = InvertedCoefficientsARMA[PacfToAR[ztaAR], 
          PacfToAR[ztaMA], papprox]; 0.5*Log[DtaARMA[ztaAR, ztaMA]] + 
         0.5*n*Log[Qd[Prepend[phiApprox, -1.], \[DoubleStruckCapitalD]]/n]]];
 
FastLoglikelihoodARMA[n_, \[DoubleStruckCapitalD]_, ztaAR_, ztaMA_, ztasAR_, 
     ztasMA_, s_] := Module[{phiApprox, papprox}, 
     -Which[Length[ztasAR] == 0 && Length[ztasMA] == 0, 
       -FastLoglikelihoodARMA[n, \[DoubleStruckCapitalD], ztaAR, ztaMA], 
       Length[ztaMA] == 0 && Length[ztasMA] == 0, 0.5*Log[Dta[ztaAR]] + 
        0.5*n*Log[Qd[Prepend[First[SARMAToARMA[PacfToAR[ztaAR], {}, PacfToAR[
                ztasAR], {}, s]], -1.], \[DoubleStruckCapitalD]]/n], 
       Length[ztaAR] == 0 && Length[ztasAR] == 0, 
       papprox = First[Dimensions[\[DoubleStruckCapitalD]]] - 1; 
        phiApprox = InvertedCoefficientsARMA[Sequence @@ SARMAToARMA[{}, 
            PacfToAR[ztaMA], {}, PacfToAR[ztasMA], s], papprox]; 
        0.5*s*Log[Dta[ztaMA]*Dta[ztasMA]] + 
         0.5*n*Log[Qd[Prepend[phiApprox, -1.], \[DoubleStruckCapitalD]]/n], 
       True, papprox = First[Dimensions[\[DoubleStruckCapitalD]]] - 1; 
        phiApprox = InvertedCoefficientsARMA[Sequence @@ 
           SARMAToARMA[PacfToAR[ztaAR], PacfToAR[ztaMA], PacfToAR[ztasAR], 
            PacfToAR[ztasMA], s], papprox]; 0.5*Log[DtaARMA[ztaAR, ztaMA]] + 
         0.5*n*Log[Qd[Prepend[phiApprox, -1.], \[DoubleStruckCapitalD]]/n]]];

HRARMA[p_, q_, x_, (opts___)?OptionQ] := 
    Module[{z, phi, res, pqmax, y, X, f, beta, vars, phiIni, thetaIni}, 
     z = x - Mean[x]; r = AROrder /. {} /. Options[HRARMA]; 
      phi = ARCoefficients /. PacfBurg[z, r, LinearPredictor -> True]; 
      res = BackcastResiduals[z, phi, {}]; pqmax = Max[p, q]; 
      y = Drop[z, pqmax]; X = Table[Take[Drop[z, pqmax - i], Length[y]], 
        {i, p}]; X = Partition[Flatten[{X, Table[Take[Drop[res, pqmax - i], 
            Length[y]], {i, q}], y}], Length[y]]; 
      vars = Table[Subscript[\[GothicCapitalU], i], {i, p + q}]; 
      f = Fit[Transpose[X], Flatten[{1, vars}], vars]; 
      beta = Table[First[f[[i]]], {i, 2, p + q + 1}]; 
      phiIni = If[p > 0, beta[[Range[p]]], {}]; 
      thetaIni = If[q > 0, beta[[p + Range[q]]], {}]; {phiIni, thetaIni}];

ResidualSDARMA[n_, \[DoubleStruckCapitalD]_, ztaAR_, ztaMA_] := 
    Module[{phiApprox, nl, papprox}, Which[Length[ztaMA] == 0 && 
       Length[ztaAR] == 0, \[DoubleStruckCapitalD][[1,1]]/n, 
      Length[ztaMA] == 0, Qd[Prepend[PacfToAR[ztaAR], -1.], 
        \[DoubleStruckCapitalD]]/n, Length[ztaAR] == 0, 
      papprox = First[Dimensions[\[DoubleStruckCapitalD]]] - 1; 
       phiApprox = InvertedCoefficientsARMA[{}, PacfToAR[ztaMA], papprox]; 
       Qd[Prepend[phiApprox, -1.], \[DoubleStruckCapitalD]]/n, True, 
      papprox = First[Dimensions[\[DoubleStruckCapitalD]]] - 1; 
       phiApprox = InvertedCoefficientsARMA[PacfToAR[ztaAR], PacfToAR[ztaMA], 
         papprox]; Qd[Prepend[phiApprox, -1.], \[DoubleStruckCapitalD]]/n]];
 
ResidualSDARMA[n_, \[DoubleStruckCapitalD]_, ztaAR_, ztaMA_, ztasAR_, 
     ztasMA_, s_] := Module[{phiApprox, papprox}, 
     Which[Length[ztasAR] == 0 && Length[ztasMA] == 0, 
      ResidualSDARMA[n, \[DoubleStruckCapitalD], ztaAR, ztaMA], 
      Length[ztaMA] == 0 && Length[ztasMA] == 0, 
      Qd[Prepend[First[SARMAToARMA[PacfToAR[ztaAR], {}, PacfToAR[ztasAR], {}, 
           s]], -1.], \[DoubleStruckCapitalD]]/n, Length[ztaAR] == 0 && 
       Length[ztasAR] == 0, 
      papprox = First[Dimensions[\[DoubleStruckCapitalD]]] - 1; 
       phiApprox = InvertedCoefficientsARMA[Sequence @@ SARMAToARMA[{}, 
           PacfToAR[ztaMA], {}, PacfToAR[ztasMA], s], papprox]; 
       Qd[Prepend[phiApprox, -1.], \[DoubleStruckCapitalD]]/n, True, 
      papprox = First[Dimensions[\[DoubleStruckCapitalD]]] - 1; 
       phiApprox = InvertedCoefficientsARMA[Sequence @@ 
          SARMAToARMA[PacfToAR[ztaAR], PacfToAR[ztaMA], PacfToAR[ztasAR], 
           PacfToAR[ztasMA], s], papprox]; Qd[Prepend[phiApprox, -1.], 
         \[DoubleStruckCapitalD]]/n]];
 
GetARMA\[DoubleStruckCapitalD][(p_)?NumberQ, (q_)?NumberQ, 
     (pApprox_)?NumberQ, (y_)?NumberVectorQ] := 
    Module[{}, Which[p == 0 && q == 0, {{y . y}}, q == 0, 
      ChampernowneD[y, p, MeanZero -> True], q > 0, ChampernowneD[y, pApprox, 
       MeanZero -> True]]];
 
GetARMA\[DoubleStruckCapitalD][(p_)?NumberQ, (q_)?NumberQ, (ps_)?NumberQ, 
     (qs_)?NumberQ, (s_)?NumberQ, (pApprox_)?NumberQ, (y_)?NumberVectorQ] := 
    Module[{}, \[DoubleStruckCapitalD] = If[q == 0 && qs == 0, 
       ChampernowneD[y, p + ps*s, MeanZero -> True], 
       ChampernowneD[y, pApprox, MeanZero -> True]]];


GetARMAMeanFit[p_, q_, pApprox_, z_, (opts___)?OptionQ] := 
    Module[{init, zinit}, init = InitialEstimates /. {} /. 
        Options[GetARMAMeanFit]; zinit = If[init === 0, Array[0. & , p + q], 
        init]; pApp = If[q == 0, p, pApprox]; y = z - (mz = Mean[z]); 
      ans = GetARMAFit[p, q, pApp, y, InitialEstimates -> init]; 
      logL = First[ans]; beta = Last[ans]; g = TacvfARMA[Sequence @@ beta, 
        pApp]; coefAR = ARCoefficients /. PacfBurg[y, pApp, 
         LinearPredictor -> True]; meanCor = getAR\[Mu][y, coefAR]; 
      x = y - meanCor; ans = GetARMAFit[p, q, pApp, x, InitialEstimates -> 
         Sequence @@ beta]; err = Abs[logL - First[ans]]; iter = 1; 
      While[err > 0.001 && iter++ < 5, logL = First[ans]; beta = Last[ans]; 
        g = TacvfARMA[Sequence @@ beta, pApp]; 
        coefAR = ARCoefficients /. PacfBurg[y, pApp, LinearPredictor -> 
            True]; meanCor = getAR\[Mu][y, coefAR]; x = y - meanCor; 
        ans = GetARMAFit[p, q, pApp, x, InitialEstimates -> 
           Sequence @@ beta]; err = Abs[logL - First[ans]]; ]; 
      {ans, Mu -> mz + meanCor, NumIterations -> iter, ErrorTolerance -> err}];

sumsqPhi[\[Phi]_] := Total[Flatten[{1, -\[Phi]}]^2];

getAR\[Mu][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[TacvfARMA[{}, phi, Length[phi]], 1]]], 1]];
 
(* THEORETICAL AUTOCORRELATION *)

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, {}}, 
       Message[TacvfARMA::singular]; Return[{}]]; 
      lagtrunc2 = LagTruncation /. {opts} /. Options[TacvfARMA]; 
      d = FractionalDifference /. {opts} /. Options[TacvfARMA]; 
      Gphi = If[Length[phi] == 0, {0.}, phi]; 
      Gtheta = If[Length[theta] == 0, {0.}, theta]; p = Length[Gphi]; 
      q = Length[Gtheta]; rp1 = 1 + Max[p, q]; qp1 = 1 + q; 
      Which[q == 0 && p == 0 && d == 0, Prepend[Table[0., {MaxLag}], 1], 
       q == 0 && p == 0 && d != 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 == 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[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;

(* INFORMATION MATRIX *)

InformationMatrixARMA[(phi_)?VectorQ, (theta_)?VectorQ, (opts___)?OptionQ] := 
    Module[{vvmatrix, uvmatrix, uumatrix, i, id}, 
     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 = {}]; 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, {{}}]];
 
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}, 
     i1 = InformationMatrixARMA[phi, theta]; 
      i2 = InformationMatrixARMA[phis, thetas]; 
      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]]; 
      Which[ps == 0 && qs == 0, i1, p == 0 && q == 0, i2, True, 
       Join[RJoin[i1, i12], RJoin[Transpose[i12], i2]]]];

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}]; 
      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; Which[i > q && i > j && imj <= q, theta[[imj]], 
       i > q && imj == 0, -1, i <= q && ijq > 0 && ijq <= p, phi[[ijq]], 
       i <= q && ijq == 0, -1, True, 0]];
 
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}]];
 
(* DIAGNOSTIC CHECKS AND SUMMARIES *)

CoefARMA[ans_] := Module[{}, {ParameterTable /. ans}//First];
 
SummaryARMA[ans_] := Module[{}, {Loglikelihood -> (Loglikelihood /. ans), 
      AIC -> (AIC /. ans), BIC -> (BIC /. ans), ARMAModelOrder -> 
       (ARMAModel /. ans), PortmanteauDiagnostic -> (PortmanteauDiagnostic /. 
        ans)}];

FullSummaryARMA[ans_]:=Drop[Drop[ans,-2],{3}];
   
(* SPECTRAL DENSITY FUNCTION *)

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

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

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, 
	  	SdfMAut[theta], 
	  	Length[theta] == 0, 
          	SdfARut[phi], 
	  	True, 
	  	SdfARut[phi]*SdfMAut[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*SdfFDWNut[d, n]]}]
   ];
 
SdfMAut[(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]]];
 
SdfARut[(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]
    ]];

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

SimulateGaussianARMA[(phi_)?InvertibleQ, (theta_)?NumberVectorQ, 
     n_Integer /; n > 0, (opts___)?OptionQ] := 
    Module[{p = Length[phi], q = Length[theta], sd}, 
     sd = Sqrt[InnovationVariance /. {opts} /. Options[SimulateARIMA]]; 
      z = sd*Which[p > 0 && q == 0, SimulateGaussianAR[phi, n], 
         p == 0 && q > 0, SimulateGaussianMA[theta, n], p > 0 && q > 0, 
         SimulateGaussianMixedARMA[phi, theta, n], p == 0 && q == 0, 
         RandomArray[NormalDistribution[0, 1], n]]];

SimulateGaussianMA[(\[Theta]_)?NumberVectorQ, n_Integer /; n > 0] := 
    ComputeMA[N[\[Theta]], RandomArray[NormalDistribution[0., 1.], 
      n + Length[\[Theta]]]];

SimulateGaussianMixedARMA[(phi_)?InvertibleQ, (theta_)?NumberVectorQ, 
     n_Integer /; n > 0] := Module[{p = Length[phi], q = Length[theta], sd, 
      svd}, v11 = Toeplitz[TacvfARMA[phi, theta, p - 1]]; 
      psi = Take[ImpulseCoefficientsARMA[phi, theta, q], q]; 
      v12 = Table[If[j - i < 0, 0, psi[[1 + j - i]]], {i, p}, {j, q}]; 
      v = BlockMatrix[{{v11, v12}, {Transpose[v12], IdentityMatrix[q]}}]; 
      svd = SingularValueDecomposition[v]; 
      zpa = RandomArray[NormalDistribution[0, 1], {p + q}] . 
        (Dot[Sequence @@ {svd[[1]], Sqrt[svd[[2]]], svd[[3]]}]); 
      zp = Reverse[Take[zpa, p]]; a = Join[Drop[zpa, p], 
        RandomArray[NormalDistribution[0, 1], n - p]]; 
      u = ComputeMA[N[theta], a]; FastGetSim[zp, u, phi]];

SimulateGaussianARIMA[(phi_)?NumberVectorQ, (theta_)?NumberVectorQ, 
     (sphi_)?NumberVectorQ, (stheta_)?NumberVectorQ, (s_)?IntegerQ, 
     (n_)?IntegerQ, (opts___)?OptionQ] := Module[{diff, gphi, gtheta, sd}, 
     sd = Sqrt[InnovationVariance /. {opts} /. Options[SimulateARIMA]]; 
      diff = Differencing /. {opts} /. Options[SimulateARIMA]; 
      {gphi, gtheta} = SARMAToARMA[phi, theta, sphi, stheta, s]; 
      w = SimulateGaussianARMA[gphi, gtheta, n, InnovationVariance -> sd]; 
      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]];

SimulateNonGaussianARMA[phi_, theta_, Q_, a_] := 
    SimMA[ImpulseCoefficientsARMA[phi, theta, Q], a];

SimMA[psi_, a_] := ListConvolve[psi, a];

(*=========================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]];
 
(* UTILITIES *)
BoundaryTestQ[(phi_)?NumberVectorQ, eps_:0] := 
    Module[{},  !Or @@ (True === #1 & ) /@ (Abs[#1 - 1] < eps & ) /@ 
         Abs[ARToPacf[phi]]];

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_, 0, s_] := z;
 
TSDifference[z_, d_ /; d < 0, s_] = Null;
 
TSDifference[z_, d_ /; d > 0, s_:1] := With[{w = TSDifference[z, d - 1, s]}, 
     Drop[w, s] - Drop[w, -s]];
 
TSDifference[z_, d_:1, s_:1] := With[{w = TSDifference[z, d - 1, s]}, 
     Drop[w, s] - Drop[w, -s]];
 
TSDifference[z_, d_:1, ds_:0, s_:1] := TSDifference[TSDifference[z, d, 1], 
     ds, s];

DtaARMA[ztaAR_, ztaMA_] := Module[{phi, theta}, 
     Which[Length[ztaAR] > 0 && Length[ztaMA] > 0, 
      phistar = GetPhiStar[PacfToAR[ztaAR], PacfToAR[ztaMA]]; 
       (Dta[ztaAR]*Dta[ztaMA])^2/DetAR[phistar], Length[ztaAR] == 0 && 
       Length[ztaMA] > 0, Dta[ztaMA], Length[ztaAR] > 0 && 
       Length[ztaMA] == 0, Dta[ztaAR], True, 1]];
 
GetPhiStar[phi_, theta_] := 
    -Drop[CoefficientList[Collect[Prepend[-phi, 1] . Table[B^(i - 1), 
           {i, Length[phi] + 1}]*Prepend[-theta, 1] . Table[B^(i - 1), 
           {i, Length[theta] + 1}], B], B], 1];
 
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];

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]]; Drop[ImpulseCoefficientsARMA[theta, gphi, 
          MaxLag]*Prepend[Table[-1, {MaxLag}], -1], 1]] /; MaxLag > 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;

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

RJoin[a_, b_] := Which[a == {{}} && b == {{}}, {}, a == {{}} && b != {{}}, b, 
     a != {{}} && b == {{}}, a, True, Transpose[Join[Transpose[a], 
       Transpose[b]]]];


(* COMPILED FUNCTIONS *)

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

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\[Equal]0,-b,
        LinearSolve[
          Table[With[{ij=i-j},
                Which[ij\[Equal]0,-1,(ij<0)||(ij>p),0,True,
                  phi\[LeftDoubleBracket]ij\[RightDoubleBracket]]]+
              If[j\[Equal]1,0.0,
                With[{ij=i+j-2},
                  Which[ij\[Equal]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]]]]];

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

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

ComputeMA=
    Compile[{{theta,_Real,1},{a,_Real,1}},
      With[{q=Length[theta]},Partition[a,q+1,1].Prepend[-theta,1.0]]];
 
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[];

