اموزش سري زماني در splus

مقدمه
يک سري زماني دنباله اي مرتب شده ازمشاهدات است.گرچه معمولا برحسب زمان , بويژه در فواصل زماني مساوي مرتب مي شود , ولي مرتب شدن ممکن باتوجه به
ابعادديگري چون فاصله نيزباشد.
يک سري زماني نظير علائم الکتريکي وولتاژ که به طورپيوسته در زمان ثبت ميشوند , يک سري پيوسته است.
يک سري زماني نظير نرخهاي بهره و حجم فروش را که فقط درفواصل مشخص در نظرگرفته ميشوند , گسسته ميگويند.
براي مطالعه سريهاي زماني , اهداف مختلفي وجود دارد , که اين اهداف درک و مکانيسم توليدي , پيش بيني مقادير آينده و کنترل بهينيک سيستم را شامل ميشود.طبيعت ذاتي يک سري زماني , وابسته يا همبسته بودن مشاهدات آن است و بنابراين روشها و فنون آماري که مبتني برفرض مستقل بودن است , ديگرکاربرد ندارد و روشهاي متفاوتي موردنياز هستند.
به روند تعيين مدل يا مدل شناسي آماري موجود , تحليل سريهاي زماني گفته ميشود.

تحليل سري زماني(الگوسازي)
در تحليل سري زماني مدلهاي مختلفي مورد بحث قرارميگيرند که ازجمله: ..,SARIMA,ARIMA,IMA ,ARI ,MA , AR که در همه اين مدلها توجه به چند نکته بسيار ضروري است.زيرا اين نکات اساس تحليل سريهاي زماني را تشکيل مي دهند .
و درغير اين صورت الگوي بدست آمده, الگوي مناسبي نخواهد بود:
1-ايستايي مدل
2-تثبيت واريانس

بررسي داده ها و هدف از تحليل آنها
داده ها مربوط به ميانگين نمونه اي وزن ترد در هرلاستيک درزمانهاي مختلف است که ازخط توليد نمونه گيري شده اند.هدف از تحليل اين داده ها بدست آوردن الگويي است که بتوان رفتار آينده سيستم را به طور کامل از دانش کنوني و محصول آينده بيان نمود.

تحليل داده ها :
قدم اول در هر تجزيه و تحليل سري زماني بررسي وضعيت ايستايي مدل است.اين بررسي معمولا به وسيله ترسيم گرافهاي pacf,acf, tsplot صورت ميگيرد.
اين گرافها در تشخيص نوع سري ازنظر ايستايي و ناايستايي , ميانگين و واريانس و نيز تشخيص اينکه داده ها داراي يک روند فصلي است يا نه, بسيار مناسبند.
دستورترسيم اين گرافها در نرم افزار بصورت زيراست:

بردار مشاهداتx=

Tsplot(x)

Acf(x)

Acf(x,type=“p”,lag=)

بنظر ميرسد که متوسط وزن ترد در لاستيکهاي نمونه گيري شده در پايان خط توليد کارخانه در حول سطح ثابتي تغيير نمي کند و شاهد روند ضعيفي در مشاهدات هستيم . يعني با افزايش دادهها در طول زمان تمرکز داده ها در حول ميانگين کاهش مي يابد و روند صعودي و گاها نزولي را مشاهده مي کنيم .اين روند را از گراف مربوط به acf در شکل يک ميتوان تشخيص داد .بدين ترتيب که acf دادههاي داراي روند به صورت خيلي ضعيفي به صفر ميل ميکند .در مدلهاي ايستا معمولا پس از مرتبه اول و دوم به صورت خيلي سريع به صفر نزديک ميشود

x_c

tsplot(x)
acf(x)
acf(x,type="p")
بنظر ميرسد که متوسط وزن ترد در لاستيکهاي نمونه گيري شده در پايان خط توليد کارخانه در حول سطح ثابتي تغيير نمي کند و شاهد روند ضعيفي در مشاهدات هستيم . يعني با افزايش دادهها در طول زمان تمرکز داده ها در حول ميانگين کاهش مي يابد و روند صعودي و گاها نزولي را مشاهده مي کنيم .اين روند را از گراف مربوط به acf در شکل يک ميتوان تشخيص داد .بدين ترتيب که acf دادههاي داراي روند به صورت خيلي ضعيفي به صفر ميل ميکند .در مدلهاي ايستا معمولا پس از مرتبه اول و دوم به صورت خيلي سريع به صفر نزديک ميشود .
با ملاحظه شکل مربوطه دادهها روند نزولي تقريبا ضعيفي را نشان مي دهند .اما از آنجايي که حذف روند اساس مدلسازي است لذا حذف روند در مشاهدات ضروري به نظر مي رسد.حذف رونددرمشاهدات راتفاضلگيري ميگويندودرsplus بادستورX1<-diff(x,lag =) صورت ميگيرد .
شکل دو گرافهاي pacf,acf, tsplot مربوط به تفاضل مرتبه اول را نشان مي دهد .

تعيين پارامتر هاي مدل:
پس از حذف روند نوبت به تعيين پارامتر هاي مدل است .اين تشخيص به وسيلهpacf,acf, ويا گرافهاي مربوط به آنها صورت ميگيرد.از آنجاييکه گرافهاي مربوطه ساختار تناوبي را نشان نميدهد لذا مدل, فصلي به نظر نميرسد و درنتيجه مدل کلي براي اين داده ها ARIMA(p,d,q) که در آن p مرتبه اتو رگرسيو و d تعداد دفعات تفاضلگيري تا حذف روند است و q مرتبه ma را نشان ميدهد .همانگونه که قبلا اشاره شد مرتبه p,q را به ترتيب ميتوان بوسيله گراف acf و pacf تشخيص داد .بدين ترتيب که فرض مربوط به شماره خطي که از باند اطمينان خارج شده است را معني دار ميدانيم (شماره خط را ميپذيريم )وآن را مرتبه qيا p ميگيريم .با توجه به شکل (ب) ميبينيم که درنقاط 1و 2 اين خطوط خارج از باند هستند و با توجه به شکل (ج)در نقاط 0و1 خارج خط هستند .لذا مرتبه p ممکن است 1يا 2و مرتبه q0 يا 1باشد .
اما در روند الگو سازي در سري زماني بايد کليه زير مدلها با استفاده از ملاکهاي موجود مورد بررسي و آزمون قرار گيرند و پس از مشاهده پيش بيني ها بهترين مدل انتخاب شود ,زير مدلها کليه ترکيبات (p,d,q) ها است که به صورت زير در نظر گرفته ميشود :

p=1,2

d=1

q=0,1

(p,d,q)=(1,0,1),(1,1,1),(2,1,0),(2,1,1)

هم دیتاهای اصلی سری زمانی.متاسفانه شکلها را کامل نشد بیاوریم

x<-c(21.4,20.35,21,20.9,20.75,20.2,20.5,21.5,21.5,22. 6,21.5,21.35,21.4,21.4,

21.7,21.2,21.3,21.2,21.8,21.2,21.5,21.35,22.7,22.5 ,21.9,22.5,21.35,21.5,2

1.2,20.85,20.9,20.85,21.5,21.1,20.9,21.1,20.65,20. 75,21.5,21.5,20.9,21.25

,21.3,20.9,21.5,19.82,20.17,20.51,20.43,20.53,20.5 6,20.89,21.41,21.81,20.

7,21.18,21.2,20,21.5,20.8,21.5,22,21.1,21.3,21.7,2 1.4,20.5,21.1,21.7,21.7

5,20.67,20.88,21.11,21.22,21.1,21.46,21.42,21.18,2 0.22,21.18,21.18,21.56,

21.28,21.58,21.6,21,21.18,20.87,19.86,21.48,20.64, 20.94,21.64,21.4,21.22,

20.46,21,21.5,20.88,20.99)

پس از برسييهاي انجام شده روي هر کدام از اين مد لها و با استفاده از ملاکهاي موجود وبهترين مدل انتخاب شد و با انجام پيش بيني هاو بررسي مانده ها مدل چندان مناسب به نظر نرسيد لذا براي بار دوم تفاضل گيري با اين هدف که شايد مدرل بهتري بدست آيد و پيش بيني ها منا سب تر باشد .

x2<-diff(x,lag=2)
شکل3گرافهاي مربوط به داده ها را پس از دو بار تفا ضل گيري نشان ميدهد .
به ترتيب قبل براي pمقادير 2و4و6 وبراي q مقادير 0و2 بدست ميايد لذا کليه زير مدلها عبارتند از :

(p,d,q)=(2,2,2),(4,2,2),(6,2,2),(2,2,0),(4,2,0),(6 ,2,0),(2,1,2),(4,1,2),(6,1,2),(2,1,0),(4,1,0),(6,1 ,0)
آنچه که از ظاهر گرا فها بر مي آيد حذف روند مشاهدات اصلي است .زيرا نقاط حول يک خط به طور تقريبا مامرکز پراکنده اند و acf ان پس از يک مرحله مشخص سريعا به صفر ميل کرده است .اين نشاندهنده حذف روند است ولذا ميتوان عمليات بعدي را براي شناسا يي الگو انجام داد .
اولين قدم براي شناسا يي مدل مناسب از بين زير مدلها داشتن يک ملاک مناسب و قوي براي تشخيص است .بهترين و مناسب ترين ملاکي که معمولا مورد استفاده قرار ميگيرد ملاک آکاييکAIC)) است که هر چه از نظر مقدار کمتر باشد نمايانگر برتري مدل است . ضمنا مدل مورد نظر بايد همگرا converg))باشد .
model1 <- list(order = c(2,2,2))
model2 <- list(order = c(4,2,2))
model3 <- list(order = c(6,2,2))
model4 <- list(order = c(2,2,0))
model5 <- list(order = c(4,2,0))
model6 <- list(order = c(6,2,0))
model7 <- list(order = c(2,1,2))
model8 <- list(order = c(4,1,2))
model9 <- list(order = c(6,1,2))
model10 <- list(order = c(2,1,0))
model11 <- list(order = c(4,1,0))
mode112 <- list(order = c(6,1,0))


براي بدست آوردن اين ملاک از دستور زير در splus استفاده ميشود:

arima.mle(x,model= model,max.iter=,max.fcal=)
model<-list(order=c(p,d,q))
که در آن به جاي (p,d,q)کليه زيرمدلهاي فوق را قرارميدهيم.
همانطور که در ضميمه مشاهده ميشود هر يک از مدلها به همراه ملاک آکاييک و همگرايي آن مشخص شده است. بنابراين با توجه به مطالب گفته شده , از بين اين مدلها مدلهاي 7و8 و 9 داراي ملاک آکاييک بيشتري نسبت به ديگر مدلها هستند ,ضمن اينکه همگرا نيز ميباشند.
براي تشخيص بهترين مدل از بين اين 3 مدل ملاکهاي ديگري بکار ميبريم.اين ملاک بررسي وضعيت مانده هاست.زيرا در مدل سري زماني , مانده ها بايد مدل نوفه سفيد يا white noise باشند.
از طرفي توزيع نرمال داشته باشند.
اين موضوع را بوسيله دستور زير انجام ميدهيم:
Arima.diag(arima.mle(x,model=model))
qqnorm(Arima.diag(arima.mle(x,model=model))$std.resid)
hist(Arima.diag(arima.mle(x,model=model))$std.resid)
ks.gof(Arima.diag(arima.mle(x,model=model))$std.resid)

که در آن به جاي (p,d,q)کليه زيرمدلهاي فوق را قرارميدهيم.
همانطور که در ضميمه مشاهده ميشود هر يک از مدلها به همراه ملاک آکاييک و همگرايي آن مشخص شده است. بنابراين با توجه به مطالب گفته شده , از بين اين مدلها مدلهاي 7و8 و 9 داراي ملاک آکاييک بيشتري نسبت به ديگر مدلها هستند ,ضمن اينکه همگرا نيز ميباشند.
براي تشخيص بهترين مدل از بين اين 3 مدل ملاکهاي ديگري بکار ميبريم.اين ملاک بررسي وضعيت مانده هاست.زيرا در مدل سري زماني , مانده ها بايد مدل نوفه سفيد يا white noise باشند.
از طرفي توزيع نرمال داشته باشند.
اين موضوع را بوسيله دستور زير انجام ميدهيم:
Arima.diag(arima.mle(x,model=model))
qqnorm(Arima.diag(arima.mle(x,model=model))$std.resid)
hist(Arima.diag(arima.mle(x,model=model))$std.resid)
ks.gof(Arima.diag(arima.mle(x,model=model))$std.resid)


که در آن به جاي model هر يک از مدلهاي 7 و 8 و 9 را قرارميدهيم.
arima7<-arima.diag(mle7)
qqnorm(arima7$std.resid)
hist(arima7$std.resid)
arima8<-arima.diag(mle8)
hist(arima8$std.resid)
qqnorm(arima8$std.resid)
arima9<-arima.diag(mle9)
hist(arima9$std.resid)
qqnorm(arima9$std.resid)



ks.gof(arima7$std.resid)
ks.gof(arima8$std.resid)
ks.gof(arima9$std.resid)

در اولين نمودار مربوط به هرکدام از ازاين مدلها(در نمودار Arima.diag ) plot مربوط به مانده ها را مشاهده ميکنيم.از انجايي که در تمام اين نمودارها روندي مشاهده نميشود , ميتوان ثابت بودن واريانس را نتيجه گرفت. لذا در اين 3 مدل واريانس ثابت است و نيازي به تثبيت واريانس ديده نميشود. Plotهاي دوم و سوم هر مدل نيزنشان ميدهد که مانده ها تقريبا از مدل نوفه سفيد پيروي ميکنند.
اما Plot چهارمي که مقدار p-value را براي ljung-Box , Chi-Square statistic را محاسبه ميکند.اين ملاک به گونه اي است که نقاط برروي آن نبايد داراي يک روند نزولي و يا صعودي باشند و يا بايد اين روند حداقل باشد.
لذا با توجه به ملاکهاي موجود مدل آخر مناسب به نظر مي رسد.چراکه نقاط روند نزولي يا صودي کمي دارند. از طرفي با مقايسه هيستو گرام مانده هاي مربوط به هرکدام از اين مدلها و آزمون نرمال بودن مانده هاي اين مدل نرمال تر هستند.

پيش بيني براي مدل انتخابي :
پس از انجام پيش بيني هايبراي مدل اول که به دو صورت انجام شد (يکبار با حذف 5 مقدار آخرين سري و بار ديگر پيش بيني مقدار براي 5 زمان بعد ) پيش بيني ها تقريبا منا سب به نظر مي رسد زيرا با محاسبه واريانس خطاي پيش بيني, مقدار واريانس کم بدست امد .

x1<-arima.forecast(x[1:100],mle9$model,n=5)
x1
arima.forecast(x[1:90],mle9$model,n=5)