Print[StyleForm[StringJoin[
        "Version 1.0. Last revision: July 29, 2003. Requires Mathematica V.5."
        ]
  ,FontSize->9,FontFamily->"Times", FontColor->
  RGBColor[0., 0.25,0.25]
   ]];

(* IAPower, Version 1.0*)
(* Author: A.I. McLeod *)
(* Copyright (c), A.I. McLeod, 2003 *)
(* Requires: Mathematica V.5. Not compatible with Version 4 or before
             due to change in how FindMinimum works. *)
(* Summary:
IAPower - Power computations for intervention analysis.
*)

BeginPackage["IAPower`"]

<<Statistics`NormalDistribution`;
<<Statistics`DescriptiveStatistics`;
<<Statistics`DataManipulation`;
<<Utilities`FilterOptions`;

\[Mu]hat;\[Phi]hat;Loglikelihood;\[Omega]hat;\[Beta]hat;\[Theta]hat;
PowerZTest;SDPowerZ;PowerLRTest;SDPowerLR;

Options[IASimpleStepFit] = {TimeSeriesModel -> "AR1"};


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

Options[MA1Fit] = {MeanZero -> False};
 
MA1Fit /: 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.";
 
SimulateAR1 /: SimulateAR1::usage = "SimulateAR1[\[Phi], n] simulates a \
series of length n from a stationary AR(1) with NID(0,1) innovations.";

SimulateMA1 /: SimulateMA1::usage = "SimulateMA1[theta,n] simulates Gaussian \
MA(1) series of length n. Takes an optional parameter, d, with default value \
0. If d is not zero, then an IMA(1) is simulated with initial starting value \
0.";
 
PulseIntervention /: PulseIntervention::usage = "PulseIntervention[T,n] generates a pulse intervention\
starting at t=T and series length n";
  
RampIntervention /: RampIntervention::usage = "RampIntervention[T,n] generates a ramp intervention\
starting at t=T and series length n";
  
StepIntervention /: StepIntervention::usage = "StepIntervention[T,n] generates a step intervention\
starting at t=T and series length n";
  
TSDifference /: TSDifference::usage = "TSDifference[z,d,s] differences or \
seasonally differences a time series d  times. The default valued 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.";
 
TacvfARMA /: TacvfARMA::singular = "Non-invertible AR or MA model or very \
close to it: the A Matrix is singular or ill-conditioned.";
 
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.";

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

FirstDifference::usage=
"FirstDifference[z] takes the first difference of data in vector z.";

Options[Sd\[Omega]ARMA] = {Difference -> 0, SeasonalPeriod -> 0, 
     SeasonalDifference -> 0, SeasonalPhi -> {}, SeasonalTheta -> {}, 
     InterventionType -> "Step"};
 
Sd\[Omega]ARMA /: Sd\[Omega]ARMA::usage = "SD\[Omega]ARIMA[T,n,phi,theta] \
computes the sd of \[Omega] in the simple step intervention model with \
ARIMA(p,q) innovations. Series length is n and intervention occurs at t=T. \
For options, see Options[Sd\[Omega]ARMA].";

Options[InformationMatrixSimpleIA] = {Difference -> 0, SeasonalPeriod -> 0, 
     SeasonalDifference -> 0, SeasonalPhi -> {}, SeasonalTheta -> {}, 
     InterventionType -> "Step"};
 
InformationMatrixSimpleIA /: InformationMatrixSimpleIA::InvalidArgs = "Invali\
d option for InterventionType. Valid InterventionType options are \"Step\", \
\"Pulse\" and \"Ramp\". Default \"Step\" will be used.";
 
InformationMatrixSimpleIA /: InformationMatrixSimpleIA::usage = "InformationM\
atrixSimpleIA[T,n, phi, theta] computes the information matrix for the \
constant term and step size coefficient in the simple step intervention with \
ARMA innovations where phi and theta are vectors of length p and q \
respectively specifying the AR and MA coefficients. The series length is n \
and the intervention occurs at t=T. Options allow for pulse and ramp \
interventions as well as ARIMA and seasonal ARIMA Models. See \
Options[InformationMatrixSimpleIA] for a list of the options and their \
defaults.";

Options[PowerFunction\[Omega]ARMA] := {Tail -> "TwoSided", Difference -> 0, 
     SeasonalPeriod -> 0, SeasonalDifference -> 0, SeasonalPhi -> {}, 
     SeasonalTheta -> {}, InterventionType -> "Step"};
 
PowerFunction\[Omega]ARMA /: PowerFunction\[Omega]ARMA::usage = "Power\[Omega\
]FunctionARMA[T,n,\[Alpha],\[Delta],phi,theta] computes the power for a step \
intervention occuring at observation number t=T. The null hypothesis is \
H0:\[Delta]=0.  ARMA(p,q) with AR(p) coefficients, phi, and MA(q) \
coefficients, theta. The size of the test is \[Alpha] and size of the \
intervention step parameter is \[Omega]=\[Delta]*Sd[z], where Sd[z] denotes \
the standard deviation of an ARMA(p,q) with AR and MA coefficients, phi and \
theta. By default, a two-tailed test is used. \
PowerFunction\[Omega]ARIMA[T,n,\[Alpha],\[Delta],phi,theta,Tail->\"LowerTail\
\"] and \
PowerFunction\[Omega]ARMA[T,n,\[Alpha],\[Delta],phi,theta,Tail->\"UpperTail\"\
] may be used to compute the one-tailed alternative. For other options, see \
Options[PowerFunction\[Omega]ARMA].";

Options[Get\[Beta]] := {TimeSeriesModel -> "AR1", Tail -> "TwoSided", 
     TestType -> "Z"};

Get\[Beta]InvalidModel /: Get\[Beta]InvalidModel::InvalidModel = 
     "Valid models are \"AR1\" and \"IMA1\". Default AR1 will be used.";
 
Get\[Beta]InvalidTest /: Get\[Beta]InvalidTest::InvalidTest = 
     "Valid TestTypes are \"LR\" and \"Z\". Default Z will be used.";
 
Get\[Beta]InvalidTailTest /: Get\[Beta]InvalidTailTest::InvalidTailTest = 
     "The LR-test is two-sided. Two-sided test assumed.";

Get\[Beta] /: Get\[Beta]::usage = "Get\[Beta][T,n,\[Alpha],beta,\[Delta],NSIM\
] computes the empirical power of two tests of size \[Alpha] of null \
hypothesis: intervention step size is zero vs the two-sided alternative. The \
size of the step intervention is \[Delta] multiplied by the standard \
deviation of the noise. The noise is by default AR(1) with parameter beta. \
For IMA(1) with parameter beta, use the option TimeSeriesModel->\"IMA1\". The \
series length is n and the intervention occurs at t=T.  One-sided tests are \
also available for the Z-test as options, Tail->\"LowerTail\" and \
Tail->\"UpperTail\". Two tests are available the Z-test and the \
likelihood-ratio test. The default is TestType->\"Z\". The other option is \
TestType->\"LR\". ";
 
Options[Plot\[Beta]] = {Tail -> "TwoSided"};
 
Plot\[Beta] /: Plot\[Beta]::usage = 
"Plot\[Beta][em,deltas,T,n,\[Alpha],beta,d] 
plots the empirical power confidence intervals and the theoretical power. \
Default is a two-tailed test. For upper-tail test, \
Plot\[Beta][em,deltas,T,n,\[Alpha],beta,d,Tail\[Rule]\"UpperTail\"].
Use d=0 for AR(1) and d=1 for MA(1). The variable em contains the
simulation results.";

Begin["`Private`"]

IASimpleStepFit[z_, T_, (opts___)?OptionQ] := Module[{indmod}, 
     indmod = TimeSeriesModel /. {opts} /. Options[IASimpleStepFit]; 
      Append[If[StringMatchQ["AR1", indmod], SimpleIAFitAR1[z, T], 
        SimpleIAFitIMA1[z, T]], TimeSeriesModel -> indmod]];

AR1MaxLogLikelihood[(z_)?NumberVectorQ, (opts___)?OptionQ] := 
    Module[{d11, d12, d22, \[Rho] = 0., \[Rho]0, etol, abserr = 100., t, rl, 
      ll1, phiopt, indzero, ETOL = 1/100000}, 
     indzero = MeanZero /. {opts} /. Options[AR1Fit]; 
      t = If[indzero, 0, Mean[z]]; d11 = getD11[z - t]; d12 = getD12[z - t]; 
      d22 = getD22[z - t, d11]; \[Rho] = Chop[getPhi[d11, d12, d22]]; 
      If[ !indav &&  !indzero, While[abserr > ETOL, 
        t = If[indzero, 0, getmean[z, \[Rho]0 = \[Rho]]]; 
         d11 = getD11[z - t]; d12 = getD12[z - t]; d22 = getD22[z - t, d11]; 
         \[Rho] = Chop[getPhi[d11, d12, d22]]; abserr = Abs[\[Rho] - \[Rho]0]/
           (Abs[\[Rho]] + ETOL); ]]; ll1 = logLikAR1[\[Rho], d11, d12, d22, 
        Length[z]]; ll1];

SimpleIAFitAR1[z_, T_] := Module[{n = Length[z]}, 
     \[GothicX] = Flatten[{Array[0 & , T - 1], Array[1 & , n - T + 1]}]; 
      \[Omega]0 = Coefficient[Fit[Transpose[{\[GothicX], z}], 
         {1, \[ScriptX]}, \[ScriptX]], \[ScriptX]]; 
      \[Omega]eps = Abs[\[Omega]0]/5; 
      ans = FindMinimum[-AR1MaxLogLikelihood[z - \[Omega]*\[GothicX]], 
        {\[Omega], \[Omega]0 - \[Omega]eps, \[Omega]0 + \[Omega]eps}]; 
      \[Omega]0 = \[Omega] /. Last[ans]; {\[Phi]0, \[Mu]0} = 
       {\[Phi]hat, \[Mu]hat} /. AR1Fit[z - \[Omega]0*\[GothicX]]; 
      {Loglikelihood -> -First[ans], \[Omega]hat -> \[Omega]0, 
       \[Beta]hat -> \[Phi]0, \[Mu]hat -> \[Mu]0}];
 
AR1Fit[(z_)?NumberVectorQ, (opts___)?OptionQ] := 
 Module[{d11, d12, d22, \[Rho] = 0., \[Rho]0, etol, abserr = 100., t, rl, 
      ll1, phiopt, indzero, ETOL = 1/100000}, 
  indzero = MeanZero /. {opts} /. Options[AR1Fit]; 
  t = If[indzero, 0, Mean[z]]; 
  d11 = getD11[z - t]; 
  d12 = getD12[z - t]; 
  d22 = getD22[z - t, d11]; 
  \[Rho] = Chop[getPhi[d11, d12, d22]]; 
  If[ !indav &&  !indzero, 
    While[abserr > ETOL, 
     t = If[indzero, 0, getmean[z, \[Rho]0 = \[Rho]]]; 
     d11 = getD11[z - t]; 
     d12 = getD12[z - t]; 
     d22 = getD22[z - t, d11]; 
     \[Rho] = Chop[getPhi[d11, d12, d22]]; 
     abserr = Abs[\[Rho] - \[Rho]0]/(Abs[\[Rho]] + ETOL); ]]; 
     ll1 = logLikAR1[\[Rho], d11, d12, d22, Length[z]]; 
     {\[Phi]hat -> \[Rho], \[Mu]hat -> t, Loglikelihood -> ll1}];
   
MA1MaxLogLikelihood[(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], mz}, 
     z = If[MeanZero /. {opts} /. Options[MA1Fit], z, z - (mz = 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}}]]]]]; If[Abs[b] < 1, fb, 
       logLMA1F[b, z]]];
  
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], mz}, 
     z = If[MeanZero /. {opts} /. Options[MA1Fit], z, z - (mz = 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}}]]]]]; {\[Theta]hat -> b, 
       \[Mu]hat -> mz, 
Loglikelihood -> If[Abs[b] < 1, fb, logLMA1F[b, z]]}];

SimpleIAFitIMA1[z_, T_] := Module[{n = Length[z], zd, xd}, 
     \[GothicX] = Flatten[{Array[0 & , T - 1], Array[1 & , n - T + 1]}]; 
      xd = FirstDifference[\[GothicX]]; zd = FirstDifference[z]; 
      \[Omega]0 = Coefficient[Fit[Transpose[{xd, zd}], {1, \[ScriptX]}, 
         \[ScriptX]], \[ScriptX]]; \[Omega]eps = Abs[\[Omega]0]/5; 
      ans = FindMinimum[-MA1MaxLogLikelihood[zd - \[Omega]*xd], 
        {\[Omega], \[Omega]0 - \[Omega]eps, \[Omega]0 + \[Omega]eps}]; 
      \[Omega]0 = \[Omega] /. Last[ans]; {\[Theta]0, \[Mu]0} = 
       {\[Theta]hat, \[Mu]hat} /. MA1Fit[FirstDifference[z] - 
          \[Omega]0*FirstDifference[\[GothicX]]]; 
      {Loglikelihood -> -First[ans], \[Omega]hat -> \[Omega]0, 
       \[Beta]hat -> \[Theta]0, \[Mu]hat -> \[Mu]0}];

logLMA1F[\[Theta]_, z_] := Module[{n = Length[z]}, 
  Check[
   detma = DetMA[\[Theta], n];
   Lz = GetLz[\[Theta], z]; 
   u = Getu0[\[Theta], Lz, detma]; 
   S = GetSumSq[\[Theta], Lz, u, n]; 
   ll = (-(1/2))*Log[detma] - (n/2)*Log[S/n]; 
   ll,Infinity 
  ]
];
 
FirstDifference[z_] := Drop[z, 1] - Drop[z, -1];

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

SimulateAR1[\[Phi]_Real /; Abs[\[Phi]] < 1, n_Integer /; n > 0] := 
    Module[{z1, a = RandomArray[NormalDistribution[0, 1], n - 1]}, 
     z1 = Random[NormalDistribution[0., 1.]]/Sqrt[1. - \[Phi]^2]; 
      FoldList[\[Phi]*#1 + #2 & , z1, a]];
 
SimulateMA1[theta_, n_, d_:0] := Module[{w}, 
     a = RandomArray[NormalDistribution[0., 1.], n + 1]; 
      w = Drop[a, 1] - theta*Drop[a, -1]; If[d == 0, w, 
      Drop[FoldList[Plus, 0, w],1]]];
 
StepIA[omega_, n_, T_] := Table[If[i < T, 0, omega], {i, n}];

TSDifference[z_, 0, s_] = z;
 
TSDifference[z_, d_ /; d < 0, s_] = Null;
 
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];
 

TacvfARMA[phi_List, theta_List, MaxLag_Integer] := 
    Module[{q = Length[theta], p = Length[phi], A, b, c, theta2, phifun, rp1, 
       qp1, d}, 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;
 
TFCalc[omega_, delta_, x_] := Module[{v, gp1 = Length[omega], 
      h = Length[delta], \[GothicG] = Plus @@ omega/(1 - Plus @@ delta)}, 
     v = Partition[Join[Table[x[[1]], {gp1 - 1}], x], gp1, 1] . 
        Reverse[omega]; If[h > 0, 
       Drop[First /@ FoldList[Flatten[{delta . #1 + #2, Drop[#1, -1]}] & , 
          Table[\[GothicG]*First[x], {h}], v], 1], v]];

InformationMatrixSimpleIA[(T_)?PostiveNumericQ, (n_)?PostiveNumericQ, 
     (phi_)?VectorQ, (theta_)?VectorQ, (opts___)?OptionQ] := 
    Module[{gp, gt, d, s, ds, sPhi, sTheta, iaType, \[Phi], \[Theta], x, b, 
       i11, i12, i22}, d = Difference /. {opts} /. 
         Options[InformationMatrixSimpleIA]; d = If[d < 0, 0, d]; 
       s = SeasonalPeriod /. {opts} /. Options[InformationMatrixSimpleIA]; 
       s = If[s < 0, 0, s]; ds = SeasonalDifference /. {opts} /. 
         Options[InformationMatrixSimpleIA]; ds = If[ds < 0, 0, ds]; 
       sPhi = SeasonalPhi /. {opts} /. Options[InformationMatrixSimpleIA]; 
       sTheta = SeasonalTheta /. {opts} /. Options[
          InformationMatrixSimpleIA]; iaType = InterventionType /. {opts} /. 
         Options[InformationMatrixSimpleIA]; 
       iaType = If[MemberQ[{"Step", "Pulse", "Ramp"}, iaType], iaType, 
         Message[InformationMatrixSimpleIA::InvalidArgs]; "Step"]; 
       {\[Phi], \[Theta]} = SARMAToARMA[phi, theta, sPhi, sTheta, s]; 
       gp = Prepend[\[Phi], -1]; gt = Prepend[\[Theta], -1]; 
       x = TSDifference[Which[iaType == "Step", StepIntervention[T, n], 
          iaType == "Pulse", PulseIntervention[T, n], iaType == "Ramp", 
          RampIntervention[T, n]], d, ds, d]; \[Kappa] = 
        Plus @@ gp/Plus @@ gt; b = TFCalc[-gp, \[Theta], x]; 
       i11 = n*\[Kappa]^2; i12 = \[Kappa]*Plus @@ b; i22 = Plus @@ (b^2); 
       {{i11, i12}, {i12, i22}}] /; n > T;

PostiveNumericQ[x_] := NumericQ[x] && x > 0;
 
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];
 
StepIntervention[T_, n_] := Flatten[{Array[0 & , T - 1], 
      Array[1 & , n - T + 1]}];
 
PulseIntervention[T_, n_] := Table[If[i == T, 1, 0], {i, n}];
 
RampIntervention[T_, n_] := Table[If[i >= T, i - T + 1, 0], {i, n}];
 
PowerFunction\[Omega]ARMA[(T_)?NumberQ, (n_)?NumberQ, (\[Alpha]_)?NumberQ, 
     (\[Delta]_)?NumberQ, (phi_)?VectorQ, (theta_)?VectorQ, 
     (opts___)?OptionQ] := Module[{Z, ty, \[Omega], sd}, 
      ty = Tail /. {opts} /. Options[PowerFunction\[Omega]ARMA]; 
       Z = Which[MemberQ[{"UpperTail", "LowerTail"}, ty], 
         Quantile[NormalDistribution[0, 1], 1 - \[Alpha]], True, 
         Quantile[NormalDistribution[0, 1], 1 - \[Alpha]/2]]; 
       \[Omega] = \[Delta]*Sqrt[First[TacvfARMA[phi, theta, 0]]]; 
       sd = Sd\[Omega]ARMA[T, n, phi, theta, FilterOptions[Sd\[Omega]ARMA, 
          {opts}]]; Which[ty == "LowerTail", CDF[NormalDistribution[0, 1], 
         -Z - \[Omega]/sd], ty == "UpperTail", 
        1 - CDF[NormalDistribution[0, 1], Z - \[Omega]/sd], True, 
        CDF[NormalDistribution[0, 1], -Z - \[Omega]/sd] + 1 - 
         CDF[NormalDistribution[0, 1], Z - \[Omega]/sd]]] /; T < n;
 
Sd\[Omega]ARMA[T_, n_, phi_, theta_, (opts___)?OptionQ] := 
    Sqrt[Last[Transpose[Inverse[InformationMatrixSimpleIA[T, n, phi, theta, 
         FilterOptions[Sd\[Omega]ARMA, {opts}]]], {1, 1}]]];
  

Get\[Beta][T_, n_, \[Alpha]_, beta_, \[Delta]_, NSIM_, (opts___)?OptionQ] := 
    Module[{tableSL = {}, tsmod, stepSize, trend, z, fullModel, west, Lfull, 
      reducedModel, Lreduced, SL, wsd, \[Beta], \[Beta]sd, ty, ta}, 
     tsmod = TimeSeriesModel /. {opts} /. Options[Get\[Beta]]; 
      ta = Tail /. {opts} /. Options[Get\[Beta]]; 
      ty = TestType /. {opts} /. Options[Get\[Beta]]; 
      If[ !MemberQ[{"AR1", "IMA1"}, tsmod], 
       Message[Get\[Beta]InvalidModel::InvalidModel]; tsmod = "AR1"]; 
      ta = If[ !MemberQ[{"TwoSided", "LowerTail", "UpperTail"}, ta], 
        Message[Get\[Beta]InvalidTail::InvalidTail]; "TwoSided", ta]; 
      ty = If[ !MemberQ[{"Z", "LR"}, ty], 
        Message[Get\[Beta]InvalidTest::InvalidTest]; "Z", ty]; 
      ta = If[ty == "LR" &&  !ta == "TwoSided", 
        Message[Get\[Beta]InvalidTailTest::InvalidTailTest]; "TwoSided", ta]; 
      stepSize = \[Delta]*Which[tsmod == "IMA1", Sqrt[1. + beta^2], True, 
         1./Sqrt[1. - beta^2]]; trend = StepIA[stepSize, n, T]; 
      Do[z = Which[tsmod == "IMA1", SimulateMA1[beta, n, 1] + trend, True, 
          SimulateAR1[beta, n] + trend]; fullModel = Which[tsmod == "IMA1", 
          IASimpleStepFit[z, T, TimeSeriesModel -> "IMA1"], True, 
          IASimpleStepFit[z, T]]; west = \[Omega]hat /. fullModel; 
        SL = If[ty == "Z", wsd = Which[tsmod == "IMA1", Sd\[Omega]ARMA[T, n, 
              {}, {beta}, Difference->1], True, Sd\[Omega]ARMA[T, n, {beta}, {}]]; 
           Z = west/wsd; Which[ta == "TwoSided", 
            2*(1 - CDF[NormalDistribution[0, 1], Abs[Z]]), ta == "UpperTail", 
            1 - CDF[NormalDistribution[0, 1], Z], ta == "LowerTail", 
            CDF[NormalDistribution[0, 1], Z]], 
          Lfull = Loglikelihood /. fullModel; reducedModel = 
            Which[tsmod == "IMA1", MA1Fit[FirstDifference[z]], True, 
             AR1Fit[z]]; Lreduced = Loglikelihood /. reducedModel; 
           1 - CDF[ChiSquareDistribution[1], 2*(Lfull - Lreduced)]]; 
        tableSL = {tableSL, SL}; , {NSIM}]; 
      \[Beta] = N[Length[Select[Flatten[tableSL], #1 <= \[Alpha] & ]]/NSIM]; 
      sd\[Beta] = Sqrt[\[Beta]*((1 - \[Beta])/NSIM)]; 
      If[ty == "Z", {PowerZTest -> \[Beta], SDPowerZ -> sd\[Beta], 
        TimeSeriesModel -> tsmod, TestType -> ty, Tail -> ta}, 
       {PowerLRTest -> \[Beta], SDPowerLR -> sd\[Beta], 
        TimeSeriesModel -> tsmod, TestType -> ty, Tail -> ta}]];
 
Plot\[Beta]Old[em_, deltas_, T_, n_, \[Alpha]_, beta_, d_, (opts___)?OptionQ] := 
    Module[{crit = 1.96, betas, SD, epi, ta}, 
     ta = Tail /. {opts} /. Options[Plot\[Beta]]; 
      ta = If[ !MemberQ[{"TwoSided", "LowerTail", "UpperTail"}, ta], 
        Message[Plot\[Beta]InvalidTail::InvalidTail]; "TwoSided", ta]; 
      {betas, SD} = Transpose[em]; {phi, theta} = Which[d == 0, {{beta}, {}}, 
        True, {{}, {beta}}]; pl = StyleForm[Which[d == 0, 
         StringForm["T=``, n=``, \!\(\[Phi]\_1\)=``", T, n, 
          NumberForm[beta, {3, 2}]], True, StringForm[
          "T=``, n=``, \!\(\[Theta]\_1\)=``", T, n, NumberForm[beta, 
           {3, 2}]]], FontWeight -> "Bold", FontSize -> 12]; 
      flab = (StyleForm[#1, FontWeight -> "Bold"] & ) /@ 
        {"\[Delta]", "\[Beta]"}; epi = {AbsoluteThickness[3], RGBColor[0, 0, 1], 
        {Line /@ Transpose[Transpose /@ {{deltas, betas - crit*SD}, 
             {deltas, betas + crit*SD}}]}}; 
      Plot[PowerFunction\[Omega]ARMA[T, n, \[Alpha], \[Delta], phi, theta, 
        Difference->d, Tail -> ta], {\[Delta], Min[deltas], Max[deltas]}, PlotStyle -> {PointSize[0.0001], 
         AbsoluteThickness[1]}, PlotRange -> {Min[0, 0.94*(betas - crit*SD)], 
         Max[1, 1.06*(betas + crit*SD)]}, AxesOrigin -> {0, 0}, 
       Frame -> True, FrameLabel -> flab, PlotLabel -> pl, 
       RotateLabel -> False, Epilog -> epi]];

Plot\[Beta][em_, deltas_, T_, n_, \[Alpha]_, beta_, d_, (opts___)?OptionQ] := 
    Module[{crit = 1.96, betas, SD, epi, ta}, 
     ta = Tail /. {opts} /. Options[Plot\[Beta]]; 
      ta = If[ !MemberQ[{"TwoSided", "LowerTail", "UpperTail"}, ta], 
        Message[Plot\[Beta]InvalidTail::InvalidTail]; "TwoSided", ta]; 
      {betas, SD} = Transpose[em]; {phi, theta} = Which[d == 0, {{beta}, {}}, 
        True, {{}, {beta}}]; pl = StyleForm[Which[d == 0, 
         StringForm[" \!\(\[Phi]\_1\)=``",  
          NumberForm[beta, {3, 2}]], True, StringForm[
          "\!\(\[Theta]\_1\)=``",  NumberForm[beta, 
           {3, 2}]]], FontWeight -> "Bold", FontSize -> 12]; 
      flab = (StyleForm[#1, FontWeight -> "Bold"] & ) /@ 
        {"\[Delta]", "\[Beta]"}; 
 epi = {AbsoluteThickness[3],  
        {Line /@ Transpose[Transpose /@ {{deltas, betas - crit*SD}, 
             {deltas, betas + crit*SD}}]}
       }; 
      Plot[PowerFunction\[Omega]ARMA[T, n, \[Alpha], \[Delta], phi, theta, 
        Difference->d, Tail -> ta], {\[Delta], Min[deltas], Max[deltas]}, PlotStyle -> {PointSize[0.0001], 
         AbsoluteThickness[1]}, PlotRange -> {Min[0, 0.94*(betas - crit*SD)], 
         Max[1, 1.06*(betas + crit*SD)]}, AxesOrigin -> {0, 0}, 
       Frame -> True, FrameLabel -> flab, PlotLabel -> pl, 
       RotateLabel -> False, Epilog -> epi]];


########Compiled Functions######################

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

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

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

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

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

logLikAR1 := Compile[{{
    phi, _Real}, {D11, _Real}, {D12, _Real}, {D22, _Real}, {n, _Integer}}, \
0.5*(Log[1 - phi^2] - n*Log[(D11 - 2*phi*D12 + phi^2*D22)/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]];

End[];
EndPackage[];
