Fungsi Transfer Intervensi ARIMA - Cara Memvisualisasikan Efek

11

Saya memiliki seri waktu bulanan dengan intervensi dan saya ingin mengukur efek intervensi ini pada hasilnya. Saya menyadari seri ini agak pendek dan efeknya belum selesai.

Data

cds <- structure(c(2580L, 2263L, 3679L, 3461L, 3645L, 3716L, 3955L, 3362L,
                   2637L, 2524L, 2084L, 2031L, 2256L, 2401L, 3253L, 2881L,
                   2555L, 2585L, 3015L, 2608L, 3676L, 5763L, 4626L, 3848L,
                   4523L, 4186L, 4070L, 4000L, 3498L),
                 .Dim=c(29L, 1L),
                 .Dimnames=list(NULL, "CD"),
                 .Tsp=c(2012, 2014.33333333333, 12), class="ts")

masukkan deskripsi gambar di sini

Metodologi

1) Seri pra-intervensi (hingga Oktober 2013) digunakan dengan auto.arimafungsi tersebut. Model yang disarankan adalah ARIMA (1,0,0) dengan mean tidak nol. Plot ACF terlihat bagus.

pre <- window(cds, start=c(2012, 01), end=c(2013, 09))

mod.pre <- auto.arima(log(pre))

# Coefficients:
#          ar1  intercept
#       0.5821     7.9652
# s.e.  0.1763     0.0810
# 
# sigma^2 estimated as 0.02709:  log likelihood=7.89
# AIC=-9.77   AICc=-8.36   BIC=-6.64

2) Mengingat alur seri penuh, respons pulsa dipilih di bawah, dengan T = Okt 2013,

masukkan deskripsi gambar di sini

yang menurut cryer dan chan bisa pas sebagai berikut dengan fungsi arimax:

mod.arimax <- arimax(log(cds), order=c(1, 0, 0),
                     seasonal=list(order=c(0, 0, 0), frequency=12),
                     include.mean=TRUE,
                     xtransf=data.frame(Oct13=1 * (seq(cds) == 22)),
                     transfer=list(c(1, 1)))
mod.arimax

# Series: log(cds) 
# ARIMA(1,0,0) with non-zero mean 
# 
# Coefficients:
#          ar1  intercept  Oct13-AR1  Oct13-MA0  Oct13-MA1
#       0.7619     8.0345    -0.4429     0.4261     0.3567
# s.e.  0.1206     0.1090     0.3993     0.1340     0.1557
# 
# sigma^2 estimated as 0.02289:  log likelihood=12.71
# AIC=-15.42   AICc=-11.61   BIC=-7.22

Sisa dari ini muncul OK:

masukkan deskripsi gambar di sini

Plot pas dan aktual:

plot(fitted(mod.arimax), col="red", type="b")
lines(window(log(cds), start=c(2012, 02)), type="b")

masukkan deskripsi gambar di sini

Pertanyaan-pertanyaan

1) Apakah metodologi ini tepat untuk analisis intervensi?

2) Dapatkah saya melihat estimasi / SE untuk komponen fungsi transfer dan mengatakan bahwa efek intervensi itu signifikan?

3) Bagaimana cara memvisualisasikan efek fungsi transfer (plot?)

4) Apakah ada cara untuk memperkirakan seberapa besar intervensi meningkatkan output setelah bulan 'x'? Saya kira untuk ini (dan mungkin # 3) Saya bertanya bagaimana cara bekerja dengan persamaan model - jika ini adalah regresi linier sederhana dengan variabel dummy (misalnya) saya bisa menjalankan skenario dengan dan tanpa intervensi dan mengukur dampaknya - tapi saya tidak yakin bagaimana cara mengerjakan model jenis ini.

MENAMBAHKAN

Per permintaan, berikut adalah residu dari dua parameter.

Pertama dari fit:

fit <- arimax(log(cds), order=c(1, 0, 0),
              xtransf=
              data.frame(Oct13a=1 * (seq_along(cds) == 22),
                         Oct13b=1 * (seq_along(cds) == 22)),
              transfer=list(c(0, 0), c(1, 0)))

plot(resid(fit), type="b")

masukkan deskripsi gambar di sini

Lalu, dari fit ini

mod.arimax <- arimax(log(cds), order=c(1, 0, 0),
                     seasonal=list(order=c(0, 0, 0), frequency=12),
                     include.mean=TRUE,
                     xtransf=data.frame(Oct13=1 * (seq(cds) == 22)),
                     transfer=list(c(1, 1))) 

mod.arimax
plot(resid(mod.arimax), type="b")

masukkan deskripsi gambar di sini

B_Miner
sumber
Apakah tidak masalah jika saya memberikan solusi menggunakan perangkat lunak SAS?
peramal
Tentu, saya akan penasaran jika Anda datang dengan model yang lebih baik.
B_Miner
Ok, modelnya sedikit lebih baik dari yang awalnya diusulkan, tetapi mirip dengan @javlacalle.
peramal

Jawaban:

12

Model AR (1) dengan intervensi yang didefinisikan dalam persamaan yang diberikan dalam pertanyaan dapat dipasang seperti yang ditunjukkan di bawah ini. Perhatikan bagaimana argumen transferdidefinisikan; Anda juga memerlukan satu variabel indikator xtransfuntuk masing-masing intervensi (denyut nadi dan perubahan sementara):

require(TSA)
cds <- structure(c(2580L, 2263L, 3679L, 3461L, 3645L, 3716L, 3955L, 3362L,
                   2637L, 2524L, 2084L, 2031L, 2256L, 2401L, 3253L, 2881L,
                   2555L, 2585L, 3015L, 2608L, 3676L, 5763L, 4626L, 3848L,
                   4523L, 4186L, 4070L, 4000L, 3498L),
                 .Dim = c(29L, 1L),
                 .Dimnames = list(NULL, "CD"),
                 .Tsp = c(2012, 2014.33333333333, 12),
                 class = "ts")

fit <- arimax(log(cds), order = c(1, 0, 0), 
              xtransf = data.frame(Oct13a = 1 * (seq_along(cds) == 22), 
                                   Oct13b = 1 * (seq_along(cds) == 22)),
              transfer = list(c(0, 0), c(1, 0)))
fit
# Coefficients:
#          ar1  intercept  Oct13a-MA0  Oct13b-AR1  Oct13b-MA0
#       0.5599     7.9643      0.1251      0.9231      0.4332
# s.e.  0.1563     0.0684      0.1911      0.1146      0.2168
# sigma^2 estimated as 0.02131:  log likelihood = 14.47,  aic = -18.94

ω0ω1coeftest

require(lmtest)
coeftest(fit)
#            Estimate Std. Error  z value  Pr(>|z|)    
# ar1        0.559855   0.156334   3.5811 0.0003421 ***
# intercept  7.964324   0.068369 116.4896 < 2.2e-16 ***
# Oct13a-MA0 0.125059   0.191067   0.6545 0.5127720    
# Oct13b-AR1 0.923112   0.114581   8.0564 7.858e-16 ***
# Oct13b-MA0 0.433213   0.216835   1.9979 0.0457281 *  
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

5%

Efek intervensi dapat dikuantifikasi sebagai berikut:

intv.effect <- 1 * (seq_along(cds) == 22)
intv.effect <- ts(
  intv.effect * 0.1251 + 
  filter(intv.effect, filter = 0.9231, method = "rec", sides = 1) * 0.4332)
intv.effect <- exp(intv.effect)
tsp(intv.effect) <- tsp(cds)

Anda dapat merencanakan efek intervensi sebagai berikut:

plot(100 * (intv.effect - 1), type = "h", main = "Total intervention effect")

Efek intervensi total

ω21ω21

Secara numerik, ini adalah perkiraan peningkatan yang dihitung pada setiap titik waktu yang disebabkan oleh intervensi pada Oktober 2013:

window(100 * (intv.effect - 1), start = c(2013, 10))
#           Jan      Feb      Mar      Apr      May Jun Jul Aug Sep      Oct
# 2013                                                              74.76989
# 2014 40.60004 36.96366 33.69046 30.73844 28.07132                         
#           Nov      Dec
# 2013 49.16560 44.64838

75%

stats::arima0.9231

xreg <- cbind(
  I1 = 1 * (seq_along(cds) == 22), 
  I2 = filter(1 * (seq_along(cds) == 22), filter = 0.9231, method = "rec", 
              sides = 1))
arima(log(cds), order = c(1, 0, 0), xreg = xreg)
# Coefficients:
#          ar1  intercept      I1      I2
#       0.5598     7.9643  0.1251  0.4332
# s.e.  0.1562     0.0671  0.1563  0.1620
# sigma^2 estimated as 0.02131:  log likelihood = 14.47,  aic = -20.94

ω20.9231xregω2

Intervensi ini setara dengan additive outlier (AO) dan perubahan sementara (TC) yang didefinisikan dalam paket tsoutliers. Anda dapat menggunakan paket ini untuk mendeteksi efek-efek ini seperti yang ditunjukkan dalam jawaban oleh @forecaster atau untuk membuat regressor yang digunakan sebelumnya. Misalnya, dalam hal ini:

require(tsoutliers)
mo <- outliers(c("AO", "TC"), c(22, 22))
oe <- outliers.effects(mo, length(cds), delta = 0.9231)
arima(log(cds), order = c(1, 0, 0), xreg = oe)
# Coefficients:
#          ar1  intercept    AO22    TC22
#       0.5598     7.9643  0.1251  0.4332
# s.e.  0.1562     0.0671  0.1563  0.1620
# sigma^2 estimated as 0.02131:  log likelihood=14.47
# AIC=-20.94   AICc=-18.33   BIC=-14.1

Edit 1

Saya telah melihat bahwa persamaan yang Anda berikan dapat ditulis ulang sebagai:

(ω0+ω1)ω0ω2B1ω2BPt

dan itu dapat ditentukan seperti yang Anda gunakan transfer=list(c(1, 1)).

Seperti ditunjukkan di bawah, parameterisasi ini mengarah, dalam hal ini, ke estimasi parameter yang melibatkan efek berbeda dibandingkan dengan parameterisasi sebelumnya. Ini mengingatkan saya pada efek outlier yang inovatif, bukannya denyut nadi plus perubahan sementara.

fit2 <- arimax(log(cds), order=c(1, 0, 0), include.mean = TRUE, 
  xtransf=data.frame(Oct13 = 1 * (seq(cds) == 22)), transfer = list(c(1, 1)))
fit2
# ARIMA(1,0,0) with non-zero mean 
# Coefficients:
#          ar1  intercept  Oct13-AR1  Oct13-MA0  Oct13-MA1
#       0.7619     8.0345    -0.4429     0.4261     0.3567
# s.e.  0.1206     0.1090     0.3993     0.1340     0.1557
# sigma^2 estimated as 0.02289:  log likelihood=12.71
# AIC=-15.42   AICc=-11.61   BIC=-7.22

Saya tidak terlalu terbiasa dengan notasi paket TSAtetapi saya berpikir bahwa efek intervensi sekarang dapat diukur sebagai berikut:

intv.effect <- 1 * (seq_along(cds) == 22)
intv.effect <- ts(intv.effect * 0.4261 + 
  filter(intv.effect, filter = -0.4429, method = "rec", sides = 1) * 0.3567)
tsp(intv.effect) <- tsp(cds)
window(100 * (exp(intv.effect) - 1), start = c(2013, 10))
#              Jan         Feb         Mar         Apr         May Jun Jul Aug
# 2014  -3.0514633   1.3820052  -0.6060551   0.2696013  -0.1191747            
#      Sep         Oct         Nov         Dec
# 2013     118.7588947 -14.6135216   7.2476455

plot(100 * (exp(intv.effect) - 1), type = "h", 
  main = "Intervention effect (parameterization 2)")

parameterisasi efek intervensi 2

Efeknya dapat digambarkan sekarang sebagai peningkatan tajam pada Oktober 2013 diikuti oleh penurunan arah yang berlawanan; maka efek intervensi lenyap dengan cepat bergantian efek positif dan negatif dari penurunan berat badan.

Efek ini agak aneh tetapi mungkin dimungkinkan dalam data nyata. Pada titik ini saya akan melihat konteks data Anda dan peristiwa yang mungkin memengaruhi data tersebut. Misalnya, apakah telah terjadi perubahan kebijakan, kampanye pemasaran, penemuan, ... yang dapat menjelaskan intervensi pada Oktober 2013. Jika demikian, apakah lebih masuk akal bahwa peristiwa ini berdampak pada data seperti yang dijelaskan sebelumnya atau seperti yang kami temukan dengan parameterisasi awal?

18.9415.42

0.9

Edit 2

ω2ω2

omegas <- seq(0.5, 1, by = 0.01)
aics <- rep(NA, length(omegas))
for (i in seq(along = omegas)) {
  tc <- filter(1 * (seq_along(cds) == 22), filter = omegas[i], method = "rec", 
               sides = 1)
  tc <- ts(tc, start = start(cds), frequency = frequency(cds))
  fit <- arima(log(cds), order = c(1, 0, 0), xreg = tc)
  aics[i] <- AIC(fit)
}
omegas[which.min(aics)]
# [1] 0.88

plot(omegas, aics, main = "AIC for different values of the TC parameter")

AIC untuk nilai omega yang berbeda

ω2=0.880.9ω2=1

ω2=0.9

ω2=0.9

tc <- filter(1 * (seq.int(length(cds) + 12) == 22), filter = 0.9, method = "rec", 
             sides = 1)
tc <- ts(tc, start = start(cds), frequency = frequency(cds))
fit <- arima(window(log(cds), end = c(2013, 10)), order = c(1, 0, 0), 
             xreg = window(tc, end = c(2013, 10)))

Prakiraan dapat diperoleh dan ditampilkan sebagai berikut:

p <- predict(fit, n.ahead = 19, newxreg = window(tc, start = c(2013, 11)))

plot(cbind(window(cds, end = c(2013, 10)), exp(p$pred)), plot.type = "single", 
     ylab = "", type = "n")
lines(window(cds, end = c(2013, 10)), type = "b")
lines(window(cds, start = c(2013, 10)), col = "gray", lty = 2, type = "b")
lines(exp(p$pred), type = "b", col = "blue")
legend("topleft",
       legend = c("observed before the intervention",
           "observed after the intervention", "forecasts"),
       lty = rep(1, 3), col = c("black", "gray", "blue"), bty = "n")

nilai yang diamati dan diperkirakan

Prakiraan pertama cocok dengan nilai-nilai yang diamati (garis titik abu-abu) yang relatif baik. Perkiraan yang tersisa menunjukkan bagaimana seri akan melanjutkan jalur ke mean asli. Interval kepercayaan masih besar, mencerminkan ketidakpastian. Karena itu kita harus berhati-hati dan merevisi model ketika data baru dicatat.

95%

lines(exp(p$pred + 1.96 * p$se), lty = 2, col = "red")
lines(exp(p$pred - 1.96 * p$se), lty = 2, col = "red")
javlacalle
sumber
Ini bagus, terima kasih! Saya punya beberapa tindak lanjut jika Anda tidak keberatan. 1) Apakah proses yang saya ikuti benar? 2) Apakah Anda menganggap kecocokan model "cukup baik" untuk menggunakan perkiraan untuk mengukur efek intervensi? 3) Haruskah saya tidak dapat menggunakan parametrization saya, yaitu transfer = daftar (c (1,1)) sebagai yang setara dan mendapatkan hasil yang cukup dekat? Contoh yang saya ikuti dari buku teks menyarankan saya harus bisa, tetapi dalam contoh ini, hasilnya tidak dekat ...
B_Miner
@ B_Miner Anda benar saya telah mengedit jawaban saya.
javlacalle
Saya setuju dengan Anda bahwa dari kedua model, parametrization pertama (mungkin dengan denyut nadi yang tidak signifikan dihapus) akan paling cocok. Mengapa dua parametrization tidak menghasilkan model yang sama, ketika saya percaya mereka harus, adalah sebuah misteri. Saya akan mengirim email kepada pengembang paket (yang juga menulis buku yang menyebutkan kesetaraan mereka).
B_Miner
Data adalah jumlah sertifikat setoran dibuka per bulan. Intervensi adalah peningkatan suku bunga rata-rata, yang melonjak mulai pada 13 Oktober. Tingkat suku bunga tetap relatif konstan sejak 13 Oktober. Tampak bagi saya bahwa setelah lonjakan, permintaan untuk produk mulai berkurang - Saya tidak yakin apakah itu akan kembali ke nilai rata-rata sebelumnya atau menetap di tingkat yang lebih tinggi (dari sebelumnya).
B_Miner
B_miner, berdasarkan data yang benar-benar tidak dapat kami simpulkan, jika permintaan akan turun ke rata-rata baru.
peramal
4

Terkadang lebih sedikit lebih banyak. Dengan 30 pengamatan di tangan saya mengirimkan data ke AUTOBOX, sebuah perangkat lunak yang telah saya bantu kembangkan. Saya menyerahkan analisis berikut ini dengan harapan mendapatkan hadiah +200 (hanya bercanda!). Saya telah merencanakan Nilai Aktual dan Bersih yang secara visual menyarankan dampak "aktivitas terbaru". masukkan deskripsi gambar di sini. Model yang dikembangkan secara otomatis ditunjukkan di sini. masukkan deskripsi gambar di sinidan di sini masukkan deskripsi gambar di sini. Residu dari seri level-shifted yang agak sederhana ini disajikan di sini masukkan deskripsi gambar di sini. Statistik model ada di sini masukkan deskripsi gambar di sini. Singkatnya ada intervensi yang dapat diidentifikasi secara empiris membuat proses ARIMA; dua pulsa dan 1 level shift masukkan deskripsi gambar di sini. Grafik Aktual / Fit dan Prakiraan lebih jauh menyoroti analisis.masukkan deskripsi gambar di sini

Saya ingin seseorang ingin melihat plot residu dari model yang ditentukan sebelumnya dan menurut saya berpotensi terlalu ditentukan.

IrishStat
sumber
Saya tidak terbiasa dengan Autobox, tetapi apakah bagian noise dari model sama dengan yang saya miliki: mean non-nol dan AR (1)?
B_Miner
Apakah keluaran ini mengatakan bahwa satu-satunya "intervensi" pada 13 Oktober hingga periode waktu saat ini adalah satu pulsa untuk 13 Oktober dan kemudian rangkaian kembali ke tingkat rata-rata normal?
B_Miner
Saya menambahkan residu dari kedua parameter. Bagi saya, sepertinya yang pertama saya daftarkan (yang semula cocok dengan javlacalle) lebih baik. Setuju?
B_Miner
1) Bagian derau adalah AR (1) dengan rerata tidak nol
IrishStat
1) Bagian derau adalah AR (1) dengan rerata tidak nol; 2) Ada 2 intervensi periode 22 dan periode 3 dan setelah 13 Oktober ia kembali ke level baru yang dimulai pada 13 September; 3) Diberi pilihan antara keduanya yang Anda sebutkan, saya setuju TETAPI saya lebih suka model AUTOBOX karena kesederhanaan dan efisiensinya. Anda dapat mengetahui lebih lanjut tentang AUTOBOX dari autobox.com/cms
IrishStat
3

R

Di bawah ini adalah kode:

cds<- structure(c(2580L, 2263L, 3679L, 3461L, 3645L, 3716L, 3955L, 
                  3362L, 2637L, 2524L, 2084L, 2031L, 2256L, 2401L, 3253L, 2881L, 
                  2555L, 2585L, 3015L, 2608L, 3676L, 5763L, 4626L, 3848L, 4523L, 
                  4186L, 4070L, 4000L, 3498L), .Dim = c(29L, 1L), .Dimnames = list(
                    NULL, "CD"), .Tsp = c(2012, 2014.33333333333, 12), class = "ts")
arimatr <- tsoutliers::tso(cds,args.tsmethod=list(d=0,D=0))
plot(arimatr)
arimatr

Di bawah ini adalah perkiraan, ada peningkatan ~ 2356,3 unit pada Oktober 2013 dengan kesalahan standar ~ 481,8 dan memiliki efek peluruhan sesudahnya. Fungsi ini secara otomatis mengidentifikasi AR (1). Saya harus melakukan beberapa iterasi dan membuat perbedaan musiman dan non musiman menjadi 0, yang tercermin dalam args.tsmethod dalam fungsi tso.

Series: cds 
ARIMA(1,0,0) with non-zero mean 

Coefficients:
         ar1  intercept       TC22
      0.5969  3034.6560  2356.2914
s.e.  0.1495   206.5202   481.7981

sigma^2 estimated as 209494:  log likelihood=-219.03
AIC=446.06   AICc=447.73   BIC=451.53

Outliers:
  type ind    time coefhat tstat
1   TC  22 2013:10    2356 4.891

Di bawah ini adalah plotnya, tsoutlier adalah satu-satunya paket yang saya tahu yang dapat mencetak perubahan sementara ini dengan baik dalam plot.

masukkan deskripsi gambar di sini

Analisis ini semoga memberikan jawaban untuk 2, 3 dan 4 pertanyaan Anda meskipun menggunakan metdeologi yang berbeda. Khususnya plot dan koefisien memberikan efek dari intervensi ini dan apa yang akan terjadi jika Anda tidak memiliki intervensi ini.

Juga berharap orang lain dapat meniru plot / analisis ini menggunakan pemodelan fungsi transfer dalam R. Saya tidak yakin apakah ini bisa dilakukan dalam R, mungkin orang lain bisa memeriksa saya tentang hal ini.

peramal cuaca
sumber
Terima kasih. Ya, plot ini adalah apa yang ingin saya lakukan dari model arimax - lihat dengan dan tanpa intervensi (dan kurangi). Saya pikir fungsi filter dalam R dapat digunakan untuk menghasilkan nilai fungsi transfer untuk setiap bulan (dan kemudian hanya memetakannya untuk memvisualisasikan) tapi saya tidak tahu bagaimana melakukannya untuk fungsi intervensi pulsa acak.
B_Miner