Initialize

Needs["LinearAlgebra`MatrixManipulation`"];

Needs["Statistics`NormalDistribution`"];

Needs["Statistics`DescriptiveStatistics`"];

<< d:/math/2003/pia.m

link = Install[
        LinkLaunch["d:/MathLink/earmaMDec18_without_reordering/Debug/earmaM.\
exe"]]

LinkPatterns[link]

SimulateFSIA[n_, T_, d_,  ω_] := 
      Module[
{z}, 
z = 
          SimulateGP[TacvfFDWN[d, n - 1]]; 
z + ω*
            Table[If[i < T, 0, 1], {i, n}]
]

\(\(SZ[n_, \ T_, NREP_, \ f_, \ ω_] := 
        Module[
{k = 0, i = 0}, 
ITMax = 
            1.25*NREP; 
fhat = 1; 
While[
            i < ITMax\  && \ i < \ NREP, 
While[
              Abs[fhat] > 0.5, 
If[
                Mod[i, 350] == 0, 
LinkClose[
                  link]; 
link\  = \ 
                  Install[
                    LinkLaunch["d:/MathLink/earmaMDec18_without_reordering/Debug/earmaM.exe"]]\ ;]; 
z = 
                SimulateFSIA[n, T, 
                  f, ω]; 
outFitted = \ 
                EarmaM[{0}, {0}, {T}, {1}, 0, 0, 0, 0, 0, 0, 1, 1, 1, z, 
                  0, {}, 0, {}, 0]; 
fhat = 
                outFitted[[1]];
]; 
sehat = 
              SEOmega[n, T, {}, {}, FractionalDifference -> f] // 
                N; 
what = 
              outFitted[[2]]; 
Zobs = what/sehat; 
k += 
              If[Abs[Zobs] > 1.96, 1, 
                0]; 
\(i++\); 
fhat = 
              1;
]; 
p = 
            k/NREP // N; 
sdp = 
            Sqrt[p*\((1 - p)\)/NREP]; 
{p, 
            sdp}
];\)\)

δs = Range[-2, 2, 0.25] // Chop;

ds = {0.0,  0.2, 0.3, 0.4};

tbl = Array[0 &, {Length[δs], Length[ds]}];

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

dsTitle = StyleForm[StringJoin["f=", ToString[#]], 
              FontColor -> RGBColor[0, 0, 1], 
              FontWeight -> "Bold"] & /@ ds;

Simulation Testing

SEOmega Testing

SEOmega[100, 50, {}, {}, FractionalDifference -> 0.4]

SEOmega[100, 50, {0.0}, {}, 
      FractionalDifference -> 0.4]

SEOmega[100, 50, {}, {}, FractionalDifference -> 0.0] // 
      N

Test out with Nile Data

LinkPatterns[link]

<< d:/math/2003/datio.m

z = ReadData["d:/MathLink/earmaM/data/nile.1"];

This is \(TraditionalForm\`z\_t = μ + ω/\((1 - δ\ B)\)\ S\  + \ a\_t/\((1 - ϕ\ B)\)\), with ϕ = 0.42, ω = -239 and δ = 0.68

EarmaM[{1}, {0}, {33}, {1}, 1, 0, 0, 0, 0, 0, 1, 0, 1, z, 
      0, {}, 0, {}, 0]

This is the simple step intervention with AR(1)

EarmaM[{0}, {0}, {33}, {1}, 1, 0, 0, 0, 0, 0, 1, 0, 1, z, 
      0, {}, 0, {}, 0]

This is the simple step intervention with fractional white noise

EarmaM[{0}, {0}, {33}, {1}, 0, 0, 0, 0, 0, 0, 1, 1, 1, z, 
      0, {}, 0, {}, 0]

Test out using regression, null case

z = SimulateFSIA[100, 50, 0.0, 0.0];

x = Table[If[i < 50, 0, 1], {i, 100}];

Fit[Transpose[{x, z}], {1, S}, S]

<<Statistics`LinearRegression`

Regress[Transpose[{x, z}], {1, S}, S]

This is the simple step intervention with AR(1)

EarmaM[{0}, {0}, {50}, {1}, 1, 0, 0, 0, 0, 0, 1, 0, 1, z, 
      0, {}, 0, {}, 0]

This is the simple step intervention with fractional white noise

EarmaM[{0}, {0}, {50}, {1}, 0, 0, 0, 0, 0, 0, 1, 1, 1, z, 
      0, {}, 0, {}, 0]

Test out using regression, big step

z = SimulateFSIA[100, 50, 0.0, 0.0];

x = Table[If[i < 50, 0, 1], {i, 100}];

z = z + 100*x;

Fit[Transpose[{x, z}], {1, S}, S]

This is the simple step intervention with AR(1)

EarmaM[{0}, {0}, {50}, {1}, 1, 0, 0, 0, 0, 0, 1, 0, 1, z, 
      0, {}, 0, {}, 0]

This is the simple step intervention with fractional white noise

EarmaM[{0}, {0}, {50}, {1}, 0, 0, 0, 0, 0, 0, 1, 1, 1, z, 
      0, {}, 0, {}, 0]

Testing SimulateFSIA

z = SimulateFSIA[100, 50, 0.0,  100]

EarmaM[{0}, {1}, {50}, {1}, 0, 0, 0, 0, 0, 0, 0, 1, 1, z, 0, {}, 0, {}, 
      0]

n = 200; T = 100; d =  .2; ω = 10.0;

z = SimulateFSIA[n, T, d, ω] + 25;

outFitted = 
      EarmaM[{0}, {1}, {T}, {1}, 0, 0, 0, 0, 0, 0, 0, 1, 1, z, 0, {}, 0, {}, 
        0]

m = Mean[Take[z, T - 1]]

Mean[Drop[z, T - 1]] - m

Testing SZ

SZ[100, 50, 5, 0.0, 0.01]

Timing Experiment

ans = Timing[SZ[100, 50, 50, 0.0, 0.0]]

First[ans]/50*1000*17/3600 Hour/Second

Plan

Try the experiment for d = 0,  0.2,  0.4

Use \(TraditionalForm\`δ = \(-3\), \ \(-2.5\), \ \(-2\), \  ... , \ \(-0.5\), \ 0, \ 0.5, \ 1, \  ... , 3\)

Range[-2, 2, 0.25] // Chop

% // Length

Main Simulation

NREP = 1000; SeedRandom[333771]; n = 50; T = 25;

\(Do[
f = 
          ds[[
            j]]; 
Do[
delta = δs[[
              i]]; 
omega = 
            delta*Sqrt[
                First[TacvfFDWN[f, 
                    1]]]; 
tbl[[i, 
              j]] = 
            SZ[n, T, NREP, f, 
              omega]; 
Print["Completed δ = ", 
            delta, ",    tbl[[i,j]]=", 
            tbl[[i, 
              j]]];
, {i, 
            Length[δs]}
]; 
Print["***Completed f = ", f];
, {j, 
          Length[ds]}
] // Timing\)

tbl

tabβ = 
        Table[First[tbl[[i, j]]], {i, 
            Length[δs]}, {j, Length[ds]}];

tabsd = 
        Table[Last[tbl[[i, j]]], {i, 
            Length[δs]}, {j, Length[ds]}];

tabβ

tabsd

TableForm[tabβ, 
      TableHeadings -> {δsTitle, dsTitle}]

crit = 1.96;

flab = (StyleForm[#1, 
                FontWeight -> "Bold"] &) /@ {"δ", "\
β"};

g = Array[0 &, Length[ds]];
\(Do[
SD = \(Transpose[tabsd]\)[[
          j]]; 
betas = \(Transpose[
            tabβ]\)[[
          j]]; 
epi = {AbsoluteThickness[
            3], {Line /@ 
              Transpose[
                Transpose /@ {{δs, betas - crit*SD}, {δs, 
                      betas + crit*SD}}]}}; 
f = 
        ds[[j]]; 
pl = 
        StyleForm[StringForm["f=``", NumberForm[f, {3, 2}]], 
          FontWeight -> "Bold", 
          FontSize -> 12]; 
g[[
          j]] = 

        Plot[omega = 
            delta*Sqrt[First[TacvfFDWN[f, 1]]]; 
PowerIA[
            100, 50, omega, {}, {}, 
            FractionalDifference -> f], {delta, \(-2\), 
            2}, 
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], 
{j, 
        Length[ds]}]\)

?Save

Save["d:/math/2003/PowerIA/Simulation/FractionalNoise/fdwn.m", {g, 
        tbl}]

Show[GraphicsArray[Transpose[{g}]]];