Simulation Experiment, IMA(1)

Introduction

In this notebook the asymptotic power is compared with the empirical power as determined by simulation for the IMA(1) case.  Both the Z-test and the likelihood-ratio test are examined.  The parameters used are discussed below.

Basic Initializations

The number of simulations,

NSIM = 1000;

Set T, the intervention point

T = 25;

Set n, the length of the time series, n > T.

n = 50;

Set the random seed so that we can repeat the results.

SeedRandom[232341];

The size of the test,

α = 0.05;

The parameters in the IMA(1) model,

thetas = {-0.75, -0.5, -0.25, 0.0, 0.25, 0.5, 
          0.75};

The parameters for the step intervention change, δ  = ω/sd(ω),

deltas = Range[-2.0, 2.0, 0.25];

Output from the simulations are saved in the directory below.  This may be need to be changed.

outDirectory = "c:/p/2003/tex/Technometrics/simulations";

The total cpu time for the Z-test and Likelihood-ratio simulations,

totTime = 0;

Suitable plot/table titles,

thetasTitle = StyleForm[StringJoin["ϕ=", ToString[#]], 
              FontColor -> RGBColor[0, 0, 1], 
              FontWeight -> "Bold"] & /@ thetas;

thetasTitle = StyleForm[StringJoin["θ=", ToString[#]], 
              FontColor -> RGBColor[0, 0, 1], 
              FontWeight -> "Bold"] & /@ thetas;

deltasTitle = StyleForm[StringJoin["δ=", ToString[#]], 
              FontColor -> RGBColor[0, 0, 1], 
              FontWeight -> "Bold"] & /@ deltas;

Basic simulation functions for AR(1) and IMA(1) IA models are provided by our notebook iapower.m

<< c:/p/2003/tex/technometrics/simulations/iapower.m

<< Statistics`ContinuousDistributions`

Main Functions

Getβ does the simulations.

Options[Getβ]

?Getβ

Plotβ produces the plots of empirical and theoretical power.

Options[Plotβ]

?Plotβ

IMA(1) , Z-test,  T=25, n=50

Testing and Timing Simulation

Testing 1

nsimTest = 100; 
    Timing[Getβ[T, n, α, 0.5, 1, nsimTest, 
        TimeSeriesModel -> "IMA1"]]

Approximate time in hours required for full simulation:

Length[thetas]*Length[deltas]*First[%]*NSIM/nsimTest/3600 
         Hour/Second // N

Testing 2

Timing[
      Getβ[T, n, α, 
        thetas[[4]], 
        deltas[[10]], 100,  
        Tail -> "TwoSided", 
        TimeSeriesModel -> "IMA1"]]

Length[thetas]*Length[deltas]*First[%]*NSIM/nsimTest/3600 
         Hour/Second // N

Main Simulation. Save Results.

Timing[
      mergedTable2SZA = 
        Table[PowerZTest /. 
            Getβ[T, n, α, 
              thetas[[i]], 
              deltas[[j]], NSIM,  
              Tail -> "TwoSided", 
              TimeSeriesModel -> "IMA1"], {j, Length[deltas]}, {i, 
            Length[thetas]}]]

87923/3600 // N

totTime = +First[%];

TableForm[mergedTable2SZA, 
      TableHeadings -> {deltasTitle, thetasTitle}]

Save Data

outFile = outDirectory <> "/ima1Z.dat";

DeleteFile[outFile]

Save[outFile, {mergedTable2SZA}];

outFileCSV = outDirectory <> "/ima1Z.csv";

Export[outFileCSV, mergedTable2SZC];

Load Data and Do Graphs, Z-test

outFile = outDirectory <> "/ima1Z.dat";

Get[outFile];

Z test

x = mergedTable2SZA;
\(\(gARZ = 
        Table[T = 25; n = 50; 

          est = Table[
              x[[i, j]], {i, 
                Length[deltas]}]; 
estSd = Sqrt[est*\((1 - est)\)/NSIM]; 

          em = Transpose[{est, estSd}]; 

          Plotβ[em, deltas, T, n, 0.05, 
            thetas[[j]], 1, \ 
            Tail -> "TwoSided"], {j, Length[thetas]}];\)\)

Likelihood-ratio Test

Main Simulation

Timing[
      mergedTable2SLRA = 
        Table[PowerLRTest /. 
            Getβ[T, n, α, 
              thetas[[i]], 
              deltas[[j]], NSIM,  
              Tail -> "TwoSided", TimeSeriesModel -> "IMA1", 
              TestType -> "LR"], {j, Length[deltas]}, {i, 
            Length[thetas]}]]

totTime += First[%];

TableForm[mergedTable2SLRA, 
      TableHeadings -> {deltasTitle, thetasTitle}]

Save AR Data

outFile = outDirectory <> "/ima1LR.dat";

DeleteFile[outFile]

Save[outFile, {mergedTable2SLRA}]

outFileCSV = outDirectory <> "/ima1LR.csv";

Export[outFileCSV, mergedTable2SZC];

Load Data and Do Graphs, Likelihood-Ratio Test

outFile = outDirectory <> "/ima1LR.dat";

Get[outFile];

Z test

x = mergedTable2SLRA;
\(\(gAR60LR = 
        Table[T = 25; n = 50; 

          est = Table[
              x[[i, j]], {i, 
                Length[deltas]}]; 
estSd = Sqrt[est*\((1 - est)\)/NSIM]; 

          em = Transpose[{est, estSd}]; 

          Plotβ[em, deltas, T, n, 0.05, 
            thetas[[j]], 1, \ 
            Tail -> "TwoSided"], {j, Length[thetas]}];\)\)