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")
Metodologi
1) Seri pra-intervensi (hingga Oktober 2013) digunakan dengan auto.arima
fungsi 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,
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:
Plot pas dan aktual:
plot(fitted(mod.arimax), col="red", type="b")
lines(window(log(cds), start=c(2012, 02)), type="b")
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")
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")
sumber
Jawaban:
Model AR (1) dengan intervensi yang didefinisikan dalam persamaan yang diberikan dalam pertanyaan dapat dipasang seperti yang ditunjukkan di bawah ini. Perhatikan bagaimana argumen
transfer
didefinisikan; Anda juga memerlukan satu variabel indikatorxtransf
untuk masing-masing intervensi (denyut nadi dan perubahan sementara):coeftest
Efek intervensi dapat dikuantifikasi sebagai berikut:
Anda dapat merencanakan efek intervensi sebagai berikut:
Secara numerik, ini adalah perkiraan peningkatan yang dihitung pada setiap titik waktu yang disebabkan oleh intervensi pada Oktober 2013:
stats::arima
xreg
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:Edit 1
Saya telah melihat bahwa persamaan yang Anda berikan dapat ditulis ulang sebagai:
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.
Saya tidak terlalu terbiasa dengan notasi paket
TSA
tetapi saya berpikir bahwa efek intervensi sekarang dapat diukur sebagai berikut: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?
Edit 2
Prakiraan dapat diperoleh dan ditampilkan sebagai berikut:
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.
sumber
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". . Model yang dikembangkan secara otomatis ditunjukkan di sini. dan di sini . Residu dari seri level-shifted yang agak sederhana ini disajikan di sini . Statistik model ada di sini . Singkatnya ada intervensi yang dapat diidentifikasi secara empiris membuat proses ARIMA; dua pulsa dan 1 level shift . Grafik Aktual / Fit dan Prakiraan lebih jauh menyoroti analisis.
Saya ingin seseorang ingin melihat plot residu dari model yang ditentukan sebelumnya dan menurut saya berpotensi terlalu ditentukan.
sumber
Di bawah ini adalah kode:
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.
Di bawah ini adalah plotnya, tsoutlier adalah satu-satunya paket yang saya tahu yang dapat mencetak perubahan sementara ini dengan baik dalam plot.
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.
sumber