पायथन समय श्रृंखला विश्लेषण

शुभ दोपहर, प्रिय पाठकों।
आज के लेख में, मैं अजगर और स्टैट्समॉडल मॉड्यूल का उपयोग करके समय श्रृंखला विश्लेषण की प्रक्रिया का वर्णन करने का प्रयास करूंगा । यह मॉड्यूल सांख्यिकीय विश्लेषण और अर्थमिति के संचालन के लिए उपकरणों और तरीकों की एक विस्तृत श्रृंखला प्रदान करता है। मैं ऐसी श्रृंखला के विश्लेषण के मुख्य चरणों को दिखाने की कोशिश करूंगा, निष्कर्ष में, हम ARIMA मॉडल का निर्माण करेंगे।
एक उदाहरण के रूप में, मास्को क्षेत्र के एक गोदाम परिसर के टर्नओवर पर वास्तविक डेटा लिया जाता है।

डाउनलोड और प्रीप्रोसेस डेटा

आरंभ करने के लिए, डेटा अपलोड करें और इसे देखें:

from pandas import read_csv, DataFrame import statsmodels.api as sm from statsmodels.iolib.table import SimpleTable from sklearn.metrics import r2_score import ml_metrics as metrics In [2]: dataset = read_csv('tovar_moving.csv',';', index_col=['date_oper'], parse_dates=['date_oper'], dayfirst=True) dataset.head() 

	Otgruzka	priemka
date_oper
2009-09-01	179,667	276,712
2009-09-02	177,670	164,999
2009-09-03	152,112	189,181
2009-09-04	142,938	254,581
2009-09-05	130,741	192,486

तो, जैसा कि आप read_csv () फ़ंक्शन को देख सकते हैं, इस मामले में, उपयोग किए गए कॉलम और इंडेक्स को निर्दिष्ट करने वाले मापदंडों को निर्दिष्ट करने के अलावा, तिथि के साथ काम करने के लिए 3 और पैरामीटर हैं। आइए हम उन पर अधिक विस्तार से ध्यान दें।
parse_dates स्तंभ नामों को DateTime प्रकार में परिवर्तित करने के लिए निर्दिष्ट करता है। यह ध्यान देने योग्य है कि यदि इस कॉलम में खाली मान हैं, तो पार्सिंग विफल हो जाएगी और टाइप ऑब्जेक्ट का एक कॉलम वापस आ जाएगा। इससे बचने के लिए, Keep_default_na = गलत पैरामीटर जोड़ें।
अंतिम डेफर्स्ट पैरामीटर पार्सिंग फ़ंक्शन को बताता है कि लाइन में पहला दिन जाने वाला पहला है, न कि इसके विपरीत। यदि यह पैरामीटर सेट नहीं है, तो फ़ंक्शन सही तरीके से तिथियों को परिवर्तित नहीं कर सकता है और स्थानों में महीने और दिन को भ्रमित कर सकता है। उदाहरण के लिए, 02/01/2013 को 02-01-2013 में बदल दिया जाएगा, जो गलत होगा।
हम शिपमेंट के मूल्यों के साथ समय श्रृंखला की एक अलग श्रृंखला का चयन करते हैं:

 otg = dataset.Otgruzka otg.head() 

date_oper
2009-09-01	179,667
2009-09-02	177,670
2009-09-03	152,112
2009-09-04	142,938
2009-09-05	130,741
नाम: Otgruzka, dtype: int64

इसलिए अब हमारे पास एक समय श्रृंखला है और हम इसके विश्लेषण के लिए आगे बढ़ सकते हैं।

समय श्रृंखला विश्लेषण

आरंभ करने के लिए, आइए हमारी श्रृंखला के ग्राफ को क्रमबद्ध करें:

 otg.plot(figsize=(12,6)) 

ग्राफ से यह देखा जा सकता है कि हमारी श्रृंखला में उत्सर्जन की एक छोटी संख्या है जो प्रसार को प्रभावित करती है। इसके अलावा, प्रत्येक दिन के लिए शिपमेंट का विश्लेषण करना पूरी तरह से सही नहीं है, क्योंकि, उदाहरण के लिए, सप्ताह के अंत या शुरुआत में ऐसे दिन होंगे जिन पर माल बाकी की तुलना में बहुत अधिक भेज दिया जाता है। इसलिए, साप्ताहिक अंतराल और उस पर लदान के औसत मूल्य पर जाने के लिए समझ में आता है, यह हमें उत्सर्जन से बचाएगा और हमारी सीमा में उतार-चढ़ाव को कम करेगा। पंडों में , इसके लिए एक सुविधाजनक पुनरावर्तन () फ़ंक्शन है, क्योंकि राउंडिंग अवधि और कुल फ़ंक्शन को पैरामीटर पास किया जाता है:

 otg = otg.resample('W', how='mean') otg.plot(figsize=(12,6)) 

जैसा कि आप देख सकते हैं, नए चार्ट में कोई उज्ज्वल उत्सर्जन नहीं है और एक स्पष्ट प्रवृत्ति है। इससे हम यह निष्कर्ष निकाल सकते हैं कि श्रृंखला स्थिर नहीं है ^[1] ।

 itog = otg.describe() otg.hist() itog 

गिनती	225
मतलब	270858.285365
एसटीडी	118371.082975
मिनट	872.857143
25%	180263.428571
50%	277898.714286
75%	355587.285714
अधिकतम	552485.142857
dtype: float64

जैसा कि विशेषताओं और हिस्टोग्राम से देखा जा सकता है, हमारी श्रृंखला कम सजातीय है और एक अपेक्षाकृत छोटा बिखराव है, जैसा कि भिन्नता के गुणांक से स्पष्ट है: जहाँ - मानक विचलन - नमूने का अंकगणित माध्य। हमारे मामले में, यह बराबर है:

 print 'V = %f' % (itog['std']/itog['mean']) 

वी = 0.437022

हम समरूपता की धारणा की पुष्टि करने के लिए वितरण नाममात्र को निर्धारित करने के लिए हर्के - बेहर परीक्षण करते हैं। ऐसा करने के लिए, एक फ़ंक्शन jarque_bera () है , जो इस आंकड़े के मूल्यों को लौटाता है:

 row = [u'JB', u'p-value', u'skew', u'kurtosis'] jb_test = sm.stats.stattools.jarque_bera(otg) a = np.vstack([jb_test]) itog = SimpleTable(a, row) print itog 

इस आंकड़े का मूल्य इंगित करता है कि वितरण की सामान्यता की अशक्त परिकल्पना को कम संभावना ( शायद> 0.05 ) के साथ खारिज कर दिया जाता है, और इसलिए हमारी श्रृंखला का सामान्य वितरण होता है।
SimpleTable () फ़ंक्शन का उपयोग आउटपुट स्वरूपित करने के लिए किया जाता है। हमारे मामले में, यह मानों की एक सरणी (2 से अधिक आयाम नहीं) और स्तंभों या पंक्तियों के नाम वाली सूची प्राप्त करता है।
कई विधियाँ और मॉडल किसी श्रृंखला की स्थिरता के बारे में मान्यताओं पर आधारित हैं, लेकिन जैसा कि पहले उल्लेख किया गया था, हमारी श्रृंखला ऐसी नहीं है। इसलिए, स्थिरता परीक्षण की जांच करने के लिए, आइए हम यूनिट की जड़ों की उपस्थिति के लिए एक सामान्यीकृत डिकी-फुलर परीक्षण करें । इसके लिए स्टैटस्मॉडल मॉड्यूल में एक एडफुलर () फ़ंक्शन है:

 test = sm.tsa.adfuller(otg) print 'adf: ', test[0] print 'p-value: ', test[1] print'Critical values: ', test[4] if test[0]> test[4]['5%']: print '  ,   ' else: print '  ,  ' 

adf: -1.38835541357
पी-मूल्य: 0.58784577297
महत्वपूर्ण मान: {'5%': -2.8753374677799957, '1%': -3.4617274344627398, '10% ': -2.5741240890815571}
एकल जड़ें हैं, पंक्ति स्थिर नहीं है


प्रस्तुत परीक्षण ने श्रृंखला के गैर-स्थिरता की धारणा की पुष्टि की। कई मामलों में, श्रृंखला के अंतर को लेने से यह अनुमति देता है। यदि, उदाहरण के लिए, श्रृंखला के पहले अंतर स्थिर हैं, तो इसे एक एकीकृत प्रथम-क्रम श्रृंखला कहा जाता है।
तो, चलो हमारी श्रृंखला के लिए एकीकृत श्रृंखला के क्रम को परिभाषित करते हैं:

 otg1diff = otg.diff(periods=1).dropna() 

ऊपर दिए गए कोड में, अंतर () फ़ंक्शन निर्दिष्ट अवधि ऑफसेट के बगल में मूल श्रृंखला में अंतर की गणना करता है। ऑफसेट अवधि को अवधि पैरामीटर के रूप में पारित किया जाता है। क्योंकि चूंकि पहला मूल्य अस्पष्ट है, तो हमें इसके लिए इसे से छुटकारा पाने की आवश्यकता है और ड्रॉपना () विधि का उपयोग किया जाता है।
स्टेशनरी के लिए परिणामी श्रृंखला देखें:

 test = sm.tsa.adfuller(otg1diff) print 'adf: ', test[0] print 'p-value: ', test[1] print'Critical values: ', test[4] if test[0]> test[4]['5%']: print '  ,   ' else: print '  ,  ' 

adf: -5.95204224907
पी-वैल्यू: 2.13583392404e-07
महत्वपूर्ण मान: {'5%': -2.8755379867788462, '1%': -3.4621857592784546, '10% ': -2.574231080806213}
एक भी जड़ें नहीं हैं, पंक्ति स्थिर है

जैसा कि आप ऊपर दिए गए कोड से देख सकते हैं, पहले अंतर की परिणामी श्रृंखला स्थिर एक के पास पहुंच गई। यह सुनिश्चित करने के लिए, हम इसे कई अंतरालों में तोड़ेंगे और सुनिश्चित करें कि चटाई बन जाएगी। विभिन्न अंतराल पर अपेक्षाएँ:

 m = otg1diff.index[len(otg1diff.index)/2+1] r1 = sm.stats.DescrStatsW(otg1diff[m:]) r2 = sm.stats.DescrStatsW(otg1diff[:m]) print 'p-value: ', sm.stats.CompareMeans(r1,r2).ttest_ind()[1] 

पी-मूल्य: 0.693072039563

उच्च पी-मूल्य हमें यह दावा करने की अनुमति देता है कि साधनों की समानता की शून्य परिकल्पना सच है, जो श्रृंखला की स्थिरता को इंगित करता है। यह सुनिश्चित करने के लिए बनी हुई है कि इसके लिए कोई प्रवृत्ति नहीं है, हम अपनी नई श्रृंखला का एक ग्राफ बनाते हैं:

 otg1diff.plot(figsize=(12,6)) 

प्रवृत्ति वास्तव में अनुपस्थित है, इसलिए पहले मतभेदों की श्रृंखला स्थिर है, और हमारी प्रारंभिक श्रृंखला पहले क्रम की एक एकीकृत श्रृंखला है ।

एक समय श्रृंखला मॉडल का निर्माण

मॉडलिंग के लिए, हम पहले कई अंतरों के लिए निर्मित ARIMA मॉडल का उपयोग करेंगे।
इसलिए, एक मॉडल बनाने के लिए हमें इसके आदेश को जानना होगा, जिसमें 2 पैरामीटर शामिल होंगे:

1. p , घटक AR का क्रम है
   
2. डी एकीकृत श्रृंखला का क्रम है
   
3. क्यू घटक एमए का क्रम है
   

पैरामीटर d है और यह 1 को तोड़ता है, यह p और q निर्धारित करने के लिए रहता है। उन्हें निर्धारित करने के लिए, हमें पहले मतभेदों के लिए ऑटोकॉरेलेशन (एसीएफ) और आंशिक रूप से ऑटोकॉरेलेशन (पीएसीएफ) कार्यों का अध्ययन करने की आवश्यकता है।
ACF हमें q निर्धारित करने में मदद करेगा, क्योंकि इसके कोरोग्राम के द्वारा, हम ऑटोकॉर्पोरेशन गुणांक की संख्या निर्धारित कर सकते हैं जो MA मॉडल में 0 से बहुत अलग हैं
PACF हमें p निर्धारित करने में मदद करेगा, क्योंकि इसके correlogram द्वारा आप AR मॉडल में 0 से भिन्न अधिकतम गुणांक संख्या निर्धारित कर सकते हैं।
संबंधित correlograms का निर्माण करने के लिए, आँकड़े फ़ंक्शन पैकेज में निम्नलिखित कार्य उपलब्ध हैं: plot_acf () और plot_pacf () । वे एसीएफ और पीएसीएफ ग्राफ प्रदर्शित करते हैं , जिसमें लैग नंबर एक्स अक्ष के साथ प्लॉट किए जाते हैं, और वाई अक्ष के साथ संबंधित कार्यों के मूल्य। यह ध्यान दिया जाना चाहिए कि कार्यों में अंतराल की संख्या महत्वपूर्ण गुणांक की संख्या निर्धारित करती है। तो, हमारे कार्य इस तरह दिखते हैं:

 ig = plt.figure(figsize=(12,8)) ax1 = fig.add_subplot(211) fig = sm.graphics.tsa.plot_acf(otg1diff.values.squeeze(), lags=25, ax=ax1) ax2 = fig.add_subplot(212) fig = sm.graphics.tsa.plot_pacf(otg1diff, lags=25, ax=ax2) 

PACF correlogram का अध्ययन करने के बाद, हम उस p = 1 का निष्कर्ष निकाल सकते हैं, क्योंकि इस पर केवल 1 अंतराल शून्य से बहुत अलग है। ACF correlogram से, हम उस q = 1 को देख सकते हैं, क्योंकि लैग 1 के बाद, फ़ंक्शन का मूल्य तेजी से गिरता है।
इसलिए, जब सभी मापदंडों को जाना जाता है, तो एक मॉडल बनाना संभव है, लेकिन इसके निर्माण के लिए हम सभी डेटा नहीं लेंगे, लेकिन केवल एक हिस्सा है। हम उस भाग से डेटा छोड़ देंगे जो हमारे मॉडल के पूर्वानुमान की सटीकता को सत्यापित करने के लिए मॉडल में नहीं आया था:

 src_data_model = otg[:'2013-05-26'] model = sm.tsa.ARIMA(src_data_model, order=(1,1,1), freq='W').fit(full_output=False, disp=0) 

मॉडल में एक स्थिर की उपस्थिति के लिए प्रवृत्ति पैरामीटर जिम्मेदार है। हम परिणामी मॉडल पर जानकारी प्राप्त करते हैं:

 print model.summary() 

जैसा कि हमारे मॉडल में इस जानकारी से देखा जा सकता है, सभी गुणांक महत्वपूर्ण हैं और हम मॉडल के मूल्यांकन के लिए आगे बढ़ सकते हैं।

विश्लेषण और मॉडल का मूल्यांकन

हम "सफेद शोर" के अनुपालन के लिए इस मॉडल के अवशेषों की जांच करेंगे, और अवशेषों के कोरलोग्राम का भी विश्लेषण करेंगे, क्योंकि इससे हमें प्रतिगमन तत्वों को निर्धारित करने में मदद मिल सकती है जो शामिल और भविष्यवाणी के लिए महत्वपूर्ण हैं।
तो पहली चीज जो हम करेंगे, वह है - लिंग - बॉक्स Q- परीक्षण की परिकल्पना का परीक्षण करने के लिए कि अवशेष यादृच्छिक हैं, अर्थात्, वे "सफेद शोर" हैं। यह परीक्षण ARIMA मॉडल के अवशेषों पर किया गया है। इस प्रकार, हमें पहले मॉडल के अवशेष प्राप्त करने और उनके लिए एसीएफ बनाने की आवश्यकता है, और फिर परिणामी गुणांक के लिए परीक्षण करें। आँकड़ेमाडल्स के साथ , यह इस तरह किया जा सकता है:

 q_test = sm.tsa.stattools.acf(model.resid, qstat=True) # resid,   , qstat=True,       - print DataFrame({'Q-stat':q_test[1], 'p-value':q_test[2]}) 

इन आंकड़ों और पी-मानों के मूल्य से संकेत मिलता है कि अवशेषों की यादृच्छिकता की परिकल्पना को अस्वीकार नहीं किया गया है, और सबसे अधिक संभावना है कि यह प्रक्रिया "सफेद शोर" का प्रतिनिधित्व करती है।
अब चलो निर्धारण के गुणांक की गणना करते हैं यह समझने के लिए कि यह मॉडल कितने प्रतिशत टिप्पणियों का वर्णन करता है:

 pred = model.predict('2013-05-26','2014-12-31', typ='levels') trn = otg['2013-05-26':] r2 = r2_score(trn, pred[1:32]) print 'R^2: %1.2f' % r2 

आर ^ 2: -0.03

हमारे मॉडल का मानक विचलन ^[2] :

 metrics.rmse(trn,pred[1:32]) 

८०,९१९.०५७३६७६४२५१२

पूर्वानुमान की औसत पूर्ण त्रुटि ^[2] :

 metrics.mae(trn,pred[1:32]) 

६३,०९२.७६३२७७६५१८९५

यह चार्ट पर हमारे पूर्वानुमान को आकर्षित करने के लिए बना हुआ है:

 otg.plot(figsize=(12,6)) pred.plot(style='r--') 

निष्कर्ष

जैसा कि आप ग्राफ से देख सकते हैं, हमारा मॉडल बहुत अच्छा पूर्वानुमान नहीं बना रहा है। यह आंशिक रूप से स्रोत डेटा में आउटलेर के कारण है, जिसे हमने पूरी तरह से हटा नहीं दिया है, साथ ही साथ आँकड़ेमाडेल पैकेज के एआरआईएमए मॉड्यूल के लिए, क्योंकि यह काफी नया है। लेख यह दिखाने के उद्देश्य से है कि वास्तव में अजगर में समय श्रृंखला का विश्लेषण कैसे किया जा सकता है। मैं यह भी नोट करना चाहूंगा कि आज जांच किए गए पैकेज में प्रतिगमन विश्लेषण के विभिन्न तरीकों को पूरी तरह से लागू किया गया है (मैं आगे के विवरण में दिखाने की कोशिश करूंगा)।
सामान्य तौर पर, छोटे अध्ययनों के लिए, सांख्यिकीमॉडल पैकेज पूरी तरह से उपयुक्त है, लेकिन गंभीर वैज्ञानिक कार्य के लिए यह अभी भी कच्चा है और कुछ परीक्षण और आंकड़े गायब हैं।

संदर्भ

1. द्वितीय Eliseev। अर्थमिति
   
2. समय श्रृंखला के मॉडल की तुलना