IAPower Documentation

Purpose and Overview of Software

The IAPower.m Mathematica package is documented in this notebook.

We review the basic results in our article, Power Computation for Intervention Analysis by McLeod and Vingilis (2003, submitted for publication) and show how these results can be extended.  The usage of all of the main functions in this package is illustrated.

The notebook Tables.nb shows how Table 1 and 2 may be constructed using the IAPower.m package. SimulationAR1.nb and SimulationAR1.nb show how the simulations were carried out and how Figures 1 and 2 were constructed.  The notebooks InformationMatrixAR1.nb and InformationMatrixIMA1.nb show how the symbolic capabilities of Mathematica were used to derive the information matrix for the Simple Step Intervention Model in the case of the AR(1) and IMA(1) models.

Intervention Analysis and the Power Computation

Intervention analysis is a widely applied time series method for quantiatively assessing the impact of an intervention on data collected over time in the form of a time series.  Intervention analysis has been widely applied in numerous areas including the social, environmental and engineering sciences.  Intervention analysis was developed in the article by Box and Tiao (1976).  The monograph by Cook and Campbell discusses the application and fundamental issues of interpretation with special emphasis on the social sciences.  Interesting applications of intervention analysis to environmetrics and stochastic hydrology are given by Hipel and McLeod (1994) and to traffic safety by Vingilis, Blefgen, Lei, Sykora, and Mann (1988).

Given a time series of length n, denoted by (Subscript[z, t]),  t = 1,  ... , n we assume that an intervention occurred and began to effect at time t = T.  It is assumed that T is known.  Before the intervention occurred it is assumed that the series is generated by an ARIMA model which may be written,

(∇^d) (Subscript[z, t]) = ξ +  + ((θ(
                  B))/(ϕ(B))) (Subscript[a, t])

where (Subscript[a, t]) are the innovations which are assumed to be normal and independently distributed with mean zero and variance (Subscript[σ, a]^2).  The AR and MA components of this model are defined by ϕ(
          B) = 1 - (Subscript[ϕ, 1]) B -  ...  (Subscript[ϕ, p]) (B^p) and θ(
          B) = 1 - (Subscript[θ, 1]) B -  ...  - (Subscript[θ, q]) (B^q) where B is the backshift operator on t and it is assumed that all roots of ϕ(B) θ(B) lie outside the unit circle.  When d = 0 this is a stationary time series.  The typical nonstationary models have d = 1 or d = 2.   The backward differencing operator is defined by (∇^d)=(((1 - B))^d).  The parameter ξ is referred to as the constant term.

The general form of the simple intervention model with ARIMA disturbances may be written,

\(TraditionalForm\`∇\^d z\_t = ξ + ω\_0\ ∇\^d ℑ\_t\%\((T)\)\  + \ \(\(θ(B)\)\/\(ϕ(B)\)\) 
          a\_t\)

where (Subscript[ℑ, t]^((T))) represents a step, pulse or ramp intervention denoted respectively by,

(Subscript[S, t]^((T)))={({{0, t < T}, {1, t ≥ T}})

(Subscript[P, t]^((T)))={({{0, t ≠  T}, {1, t = T}})

(Subscript[R, t]^((T)))={({{0, t < T}, {t - T + 1, t ≥ T.}})

Under suitable regularity conditions described in Pierce (1972), the  maximum likelihood estimates for the parameter vector (ξ, (Subscript[ω, 0])) are consistent and approximately normally distributed in large samples with mean (ξ, (Subscript[ω, 0])) and variance (𝕀^(-1)), where 𝕀 is the Fisher large-sample information matrix which may be written,

𝕀=(({{n (κ^2) , \(\(Σ\_t\) κ\ b\_t\)}, {(Subscript[Σ, t]) κ (Subscript[b, t]), \(\(Σ\_t\)\(\ \)\(b\_t\%2\)\(\ \)\)}}))

where κ = -ϕ(1)/θ(1) and (Subscript[b, t]) = -ϕ(B)/θ(
              B)  (Subscript[x, t]), where \(TraditionalForm\`x\_t = \(\(∇\^d ℑ\_t\%\((T)\)\)\(.\)\)\)

In order to detect if the intervention has caused a changed in the time series we can test the null hypothesis, (Subscript[ℍ, 0]) : (Subscript[ω, 0]) = 
        0 vs its two-sided alternative of simple negation.  After fitting the intervention analysis model (Subscript[ℍ, 0]) may be tested by computing its p-value given \(TraditionalForm\`\(\(p = 
        Pr {\(|\)\(\(ω\&^\)\_0\)\(|\)\(\(>\)\(𝒵\_\(1 - α/2\)\ \(σ\&^\)\_\(ω\_0\)\)\)}\)\(,\)\)\) where \(TraditionalForm\`\(σ\&^\)\_\(ω\_0\)\) is the estimated standard deviation of \(TraditionalForm\`\(ω\&^\)\_0\) and \(TraditionalForm\`𝒵\_\(1 - α/2\) = Φ(1 - α/2)\), where Φ(▪) denotes the cumulative distribution function of the standard normal distribution.  The power, defined as the probably of rejecting (Subscript[ℍ, 0]) when (Subscript[ω, 0]) ≠ 0 may be expressed \(TraditionalForm\`Pr {\(|\)\(\(ω\&^\)\_0\)\(|\)\(\(>\)\(𝒵\_\(1 - α/2\)\ σ\_\(ω\_0\)\)\)\(|\)\(ω\_0\)}\).  It is more convenient to work with \(TraditionalForm\`δ = ω\_0/σ\_\(ω\_0\)\) which makes the power function functionally independent of (Subscript[σ, a]^2).  The two-sided power function may then be written,

\(TraditionalForm\`β(δ) = Φ(\(-Z\_\(1 - α/2\)\)\  - δ\ σ/σ\_\(ω\_0\)) + 
        1 - Φ(
          Z\_\(1 - α/2\)\  - δ\ σ/σ\_\(ω\_0\))\)

In some cases the alternative hypothesis is not simple negation but rather (Subscript[ℍ, 1]) : (Subscript[ω, 0]) > 
          0.  For this upper-tailed test, the power function is

\(TraditionalForm\`\(β\_u\)(δ) = 
      1 - Φ(
          Z\_\(1 - α\)\  - δ\ σ/σ\_\(ω\_0\))\)

and similarly for the lower-tail test,

\(TraditionalForm\`\(β\_l\)(δ) = Φ(\(-Z\_\(1 - α\)\)\  - δ\ σ/σ\_\(ω\_0\))\).

In our article, Power Computation for Intervention Analysis, we show how these power functions may easily be computed with a hand calculator in the case of the AR(1) and IMA(1) models.  More general computations can be done with our package IAPower.m.

It is very important to understand the power computations must be made prior to collecting and fitting the model.  Power computations have sometimes been mis-applied in a mis-guided attempt to assess the usefulness of a statistical significance test after it has already been applied.  This idea is discussed in detail by Hoenig and Heisey (2001).

References

Box, G.E.P. and Tiao, G.C. (1975), Intervention analysis with applications to economic and environmental problems, Journal of the American Statistical Association 70, 70-79.

Cook, T.D. and Campbell, D.T. (1979), Quasi-Eperimentation: Design & Analysis Issues for Field Settings, Houghton Mifflin, Boston.

Hipel, K.W. and McLeod, A.I. (1994). Time Series Modelling of Water Resources and Environmental Systems.  Amesterdam: Elsevier.

Hoenig, J. M. and Heisey, D. M. (2001). The abuse of power the pervasive fallacy of power calculations for data analysis.  The American Statistician   55  19--24.

Pierce, D.A. (1972).  Least squares estimation in dynamic-disturbance time series models.  Biometrika 59, 73-78.

Vingilis, E., Blefgen, H., Lei, H., Sykora, K. and Mann, R. (1988).  An evaluation of the deterrent impact of Ontario's 12-hour licence suspension law.  Accident Analysis and Prevention 20, 9-17.

Summary of Package

First we load the package,

<< c:/p/2003/tex/Technometrics/Webpage/iapower.m

and list all the functions and variables,

?IAPower`*

Power Computation

PowerFunctionωARMA may be used to compute β(δ), (Subscript[β, u])(δ) or (Subscript[β, l])(δ) for any intervention model defined by eqns. (4)-(6) above.

?PowerFunctionωARMA

Options[PowerFunctionωARMA]

Illustrative Examples

Suppose we have T = 50, n = 100 and an ARIMA(1,0,1) model with (Subscript[ϕ, 1]) = 0.8 and (Subscript[θ, 1]) = 0.4 then for a two-sided test at level α = 0.05 the value of β(δ) with δ = 1 is given by,

PowerFunctionωARMA[50, 100, 0.05, 1.0, {0.8}, {0.4}]

The power for a corresponding upper-tail test,

PowerFunctionωARMA[50, 100, 0.05, 1.0, {0.8}, {0.4}, 
      Tail -> "UpperTail"]

For an ARIMA(1,1,1) with the same values for the other parameters the power is greatly reduced,

PowerFunctionωARMA[50, 100, 0.05, 1.0, {0.8}, {0.4}, 
      Difference -> 1]

PowerFunctionωARMA[50, 100, 0.05, 1.0, {0.8}, {0.4}, 
      Tail -> "UpperTail", Difference -> 1]

Finally lets take the multiplicative seasonal ARIMA (0, 1, 1) (Subscript[((0, 1, 1)), 12]) with (Subscript[θ, 1]) = 0.8 and (Subscript[Θ, 1]) = 0.8,

PowerFunctionωARMA[50, 100, 0.05, 1.0, {}, {0.8}, 
      Difference -> 1, SeasonalDifference -> 1, 
      SeasonalPeriod -> 12, SeasonalTheta -> {0.8}]

PowerFunctionωARMA[50, 100, 0.05, 1.0, {}, {0.8}, 
      Tail -> "UpperTail", Difference -> 1, 
      SeasonalDifference -> 1, SeasonalPeriod -> 12, 
      SeasonalTheta -> {0.8}]

Simulation of Empirical Power

The function for determining the empirical power is only for ARIMA(1,0,0) or ARIMA(0,1,1) models.

?Getβ

In the examples below we do 100 simulations exercising various options.

SetSeed[4312];

Timing[Getβ[36, 60, 0.05, 0.5, 1.0, 100]]

Timing[
      Getβ[36, 60, 0.05, 0.5, 1.0, 100, 
        TestType -> "LR"]]

Timing[
      Getβ[36, 60, 0.05, 0.5, 1.0, 100, 
        Tail -> "UpperTail"]]

Timing[
      Getβ[36, 60, 0.05, 0.5, 1.0, 100, 
        Tail -> "LowerTail"]]

Timing[
      Getβ[36, 60, 0.05, 0.5, 1.0, 100, 
        TimeSeriesModel -> "IMA1"]]

Timing[
      Getβ[36, 60, 0.05, 0.5, 1.0, 100, TestType -> "LR", 
        TimeSeriesModel -> "IMA1"]]

Timing[
      Getβ[36, 60, 0.05, 0.5, 1.0, 100, Tail -> "UpperTail", 
        TimeSeriesModel -> "IMA1"]]

Timing[
      Getβ[36, 60, 0.05, 0.5, 1.0, 100, Tail -> "LowerTail", 
        TimeSeriesModel -> "IMA1"]]

Simulated Intervention Models

Simple intervention models with ARIMA(1, 0, 0) or ARIMA(0, d, 1), where d = 0, 1 may be simulated.  First here are the two functions for simulating the disturbance component,

?SimulateAR1

?SimulateMA1

For example to simulate an ARIMA(0,1,1), with n = 20 and (Subscript[θ, 1]) = 0.8,

z = SimulateMA1[0.8, 20, 1]

Three functions are provided for various types of interventions,

?StepIntervention

StepIntervention[25, 50]

?PulseIntervention

PulseIntervention[25, 50]

?RampIntervention

RampIntervention[25, 50]

So to simulate a Simple Step Intevention with (Subscript[ω, 0]) = 2.5 and an ARIMA(0,1,1) disturbance component, with T = 12, n = 20 and (Subscript[θ, 1]) = 0.8,

z = SimulateMA1[0.8, 20, 1] + 2.5*StepIntervention[12, 20]

Fitting Models

We provide IASimpleStepFit, MA1Fit, AR1Fit and SdωARMA.  Using these functions one can fit by exact maximum likelihood the ARIMA(1,0,0) and ARIMA(0,1,1) models as well as simple step, pulse or ramp interventions with these types of disturbances.  Likelihood-ratio tests may be done or by using SdωARMA one can obtain \(TraditionalForm\`\(σ\&^\)\_\(ω\_0\)\) and compute the Z-test.

Illustrative Example

We simulate a simple step intervention model with (Subscript[ω, 0]) = 0.65 and AR(1) disturbances with ϕ = 0.5 and n = 200 and T = 100 and compute the maximum likelihood estimates.  The significance of the intervention is tested using the Z-test and a likelihood ratio test.

In order to discuss the results we will set the random seed.

SetSeed[4312];

ϕ = 0.5; n = 200; T = 100; stepSize = 0.75; 
    z = SimulateAR1[ϕ, n] + stepSize*StepIntervention[T, n];

fullModel = IASimpleStepFit[z, T]

west = ωhat /. fullModel

Lfull = Loglikelihood /. fullModel

Lreduced = Loglikelihood /. AR1Fit[z]

For the likelihood ratio test we need to load the Standard Add-On Library:

Needs["Statistics`ContinuousDistributions`"]

Now we compute the achieved p-value using a likelihood-ratio test,

SLLRTest = 
      1 - CDF[ChiSquareDistribution[1], 2*(Lfull - Lreduced)]

and for comparison the p-value using the Z-test

wsd = SdωARMA[T, n, {ϕ}, {}]

SLZtestTwoTail = 
      2*(1 - CDF[NormalDistribution[0, 1], Abs[west/wsd]])

SLZtestUpperTail = 
      1 - CDF[NormalDistribution[0, 1], west/wsd]