51R語言機(jī)器學(xué)習(xí)：時(shí)間序列與因果關(guān)系

51-R語言機(jī)器學(xué)習(xí)：時(shí)間序列與因果關(guān)系《精通機(jī)器學(xué)習(xí)：基于R第二版》學(xué)習(xí)筆記1、數(shù)據(jù)理解與數(shù)據(jù)準(zhǔn)備使用的數(shù)據(jù)記錄的是年度異常數(shù)據(jù)，它由某一年的年度地表溫度中位數(shù)與參考年度(1961~1990)平均溫度的差構(gòu)成。climate<-read.csv("./data_set/data-master/climate.csv",stringsAsFactors=F)str(climate)##'data.frame':95obs.of3variables:##$Year:int1919192019211922192319241925192619271928...##$CO2:int80693280384597096397598310621065...##$Temp:num-0.272-0.241-0.187-0.301-0.272-0.292-0.214-0.105-0.208-0.206...>summary(climate)##YearCO2########TempMin.:1919Min.:803Min.:-0.348001stQu.:1942Median:1966Mean:19661stQu.:13161stQu.:-0.13150Median:3288Median:0.00000##YearCO2########TempMin.:1919Min.:803Min.:-0.348001stQu.:1942Median:1966Mean:19661stQu.:13161stQu.:-0.13150Median:3288Median:0.00000Mean:3852Mean:0.04773####3rdQu.:1990Max.:20133rdQu.:60813rdQu.:0.18550Max.:9776Max.:0.55600>#將其轉(zhuǎn)換為時(shí)間序列結(jié)構(gòu)>climate<-ts(climate[,2:3],frequency=1,start=1919,end=2013)首先繪制兩個(gè)時(shí)間序列:library(pacman)p_load(forecast,tseries,ggfortify)plot(climate)檸掙型報(bào)靜掙型報(bào)靜:竦型報(bào)靜:陽Ji報(bào)靜掙型報(bào)靜掙型報(bào)靜:陽Ji報(bào)靜掙型報(bào)靜:陽妹2再climate變量時(shí)間序列從圖中可以看出，溫度的異常變化大約發(fā)生在20世紀(jì)70年代中期，那時(shí)開始急劇升高。再看看兩個(gè)變量之間的相關(guān)性：##C02Temp##CO21.00000000.8404215##Temp0.84042151.0000000兩個(gè)序列是高度相關(guān)的。繪制兩個(gè)序列的ACF圖和PACF圖，再看看序列結(jié)構(gòu):p_load(patchwork)>p1<-autoplot(acf(climate[,2],plot=F),main="TempACF")p2<-autoplot(pacf(climate[,2],plot=F),main="TempPACF")p3<-autoplot(acf(climate[,1],plot=F),main="CO2ACF")p4<-autoplot(pacf(climate[,1],plot=F),main="CO2PACF")>>p1+p2+p3+p4+plot_layout(ncol=2,nrow=2)ACF圖和PACF圖ACF圖中虛線表示顯著相關(guān)性的置信帶，任何一條超出置信帶上界和下界的豎線都被認(rèn)為是顯著的相關(guān)性。PACF計(jì)算的是條件相關(guān)性。ACF模式是逐漸衰減的，PACF模式則是快速衰減的?？梢约僭O(shè)這兩個(gè)序列都是自回歸的，盡管溫度中看上去存在明顯的MA項(xiàng)。下一步檢查交叉相關(guān)函數(shù)。#函數(shù)規(guī)定x要放到y(tǒng)的前面Ccf(climate[,1],climate],2],main="CCF")CCF圖CCF圖展示了溫度和二氧化碳延遲序列之間的相關(guān)性。如果x變量的負(fù)延遲序列具有強(qiáng)相關(guān)性，則x領(lǐng)先于y;如果x變量的正延遲序列具有強(qiáng)相關(guān)性，則x滯后于y。此處可以看到，二氧化碳既是個(gè)領(lǐng)先變量，也是個(gè)滯后變量。檢驗(yàn)數(shù)據(jù)是否平穩(wěn)：#擴(kuò)展迪基-福勒檢驗(yàn)adf.test(climate[,1])####AugmentedDickey-FullerTest####data:climate],1]##Dickey-Fuller=-1.1519,Lagorder=4,p-value=0.9101##alternativehypothesis:stationary>adf.test(climate[,2])####AugmentedDickey-FullerTest####data:climate],2]##Dickey-Fuller=-1.8106,Lagorder=4,p-value=0.6546##alternativehypothesis:stationary可以看到，對(duì)于這兩個(gè)序列，p值都是不顯著的，所以不能拒絕原假設(shè)(數(shù)據(jù)是不穩(wěn)定的)。2、模型構(gòu)建與模型評(píng)價(jià)主要有三項(xiàng)任務(wù)：第一，建立一個(gè)僅應(yīng)用于溫度數(shù)據(jù)的單變量預(yù)測(cè)模型；第二，基于溫度數(shù)據(jù)本身和二氧化碳排放量數(shù)據(jù)建立一個(gè)溫度數(shù)據(jù)的回歸模型；第三，使用這個(gè)模型的輸出揭示二氧化碳排放量和地表溫度異常之間是否存在格蘭杰因果關(guān)系。2.1單變量時(shí)間序列分析為地表溫度建立一個(gè)單變量預(yù)測(cè)模型，時(shí)間從二戰(zhàn)之后開始:temp<-climate],2]>>#拆分為訓(xùn)練集和測(cè)試集train<-window(temp,start=1946,end=2003)test<-window(temp,start=2004)指數(shù)平滑模型對(duì)過去的觀測(cè)進(jìn)行加權(quán)。但和移動(dòng)平均模型不同的是，在指數(shù)平滑模型中，相對(duì)于更遠(yuǎn)期的觀測(cè)來說，越是近期的觀測(cè)，所得的權(quán)重越咼?，F(xiàn)在，建立沒有阻尼趨勢(shì)的holt模型：#預(yù)測(cè)階段數(shù)量h初始狀態(tài)值的方法('optimal'或者'simple')#'optimal'表示算法會(huì)找出最優(yōu)的初始值以及平滑參數(shù)#'simple'表示使用開始的幾個(gè)觀測(cè)來計(jì)算初始值fit.holt<-holt(train,h=10,initial="optimal")>>#繪制預(yù)測(cè)值，看看模型在樣本外數(shù)據(jù)上的預(yù)測(cè)效果plot(forecast(fit.holt))lines(test,type="o")ForecastsfromHoltsmethod預(yù)測(cè)圖圖中的預(yù)測(cè)值表現(xiàn)出線性上升趨勢(shì)。下面加上阻尼趨勢(shì):fit.holtd<-holt(train,h=10,initial="optimal",damped=T)plot(forecast(fit.holtd),main="HOLTDamped")lines(test,type="o")加上阻尼趨勢(shì)的預(yù)測(cè)圖最后再建立一個(gè)ARIMA模型:##Series:train##ARIMA(0,1,1)withdrift####Coefficients:##ma1drift##-0.69490.0094##s.e.0.10410.0047####sigma八2estimatedas0.01273:loglikelihood=44.18##AIC=-82.37AICc=-81.92BIC=-76.24####Trainingseterrormeasures:##MERMSEMAEMPEMAPEMASE##Trainingset-0.00070389220.10985610.0901770726.99009153.25560.878048##ACF1##Trainingset0.1282754函數(shù)選擇的模型是MA=1、I=1,也就是帶有漂移項(xiàng)(等價(jià)于截距)的ARIMA(0,1,1)。具有一階自回歸(AR(1))和一階差分(1(1))的ARIMA模型可以表示為ARIMA(1,1,0)。通過圖形檢查它在測(cè)試數(shù)據(jù)上的表現(xiàn)：plot(forecast(fit.arima,h=10))lines(test,type="o")ForecastsfromARIMA(0T1,1)witldrift模型在測(cè)試集上的表現(xiàn)這個(gè)圖與不帶阻尼趨勢(shì)的holt方法非常相似。通過下面的代碼，可以為每種模型打分，找出具有最低誤差，即平均絕對(duì)百分誤差的模型:##[1]0.2864105##[1]0.2864105mape.holtd<-sum(abs((test-fit.holtd$mean)/test))/10print(mape.holtd)##[1]0.2211337mape.arima<-sum(abs((test-forecast(fit.arima,h=10)$mean)/test))/10print(mape.arima)##[1]0.1034813##[1]0.1034813與holt方法相比，ARIMA(0,1,1)的預(yù)測(cè)誤差要稍小一些。很顯然，不帶阻尼趨勢(shì)的模型表現(xiàn)最差。通過統(tǒng)計(jì)檢驗(yàn)和可視化證據(jù)，單變量預(yù)測(cè)模型的最好選擇似乎是ARIMA模型。2.2檢查因果關(guān)系首先演示殘差受到自相關(guān)(也稱為序列相關(guān))影響而導(dǎo)致的虛假線性回歸。然后，研究?jī)煞N實(shí)現(xiàn)格蘭杰因果關(guān)系模型的不同方法。第一種是傳統(tǒng)方法，其中兩個(gè)序列都是平穩(wěn)的。此后看看Toda和Yamamoto(1995)提出的方法，這種方法可以應(yīng)用于原始數(shù)據(jù)(有時(shí)也稱為“水平”)。2.2.1線性回歸fit.lm<-lm(Temp~CO2,data=climate)summary(fit.lm)####Call:##lm(formula=Temp~CO2,data=climate)####Residuals:##Min1QMedian3QMax##-0.36411-0.089860.000110.094750.28763####Coefficients:##EstimateStd.ErrortvaluePr(>|t|)##(Intercept)-2.430e-012.357e-02-10.31<2e-16***##CO27.548e-055.047e-0614.96<2e-16***##---##Signif.codes:0'***'0.001'**'0.01'*'0.05'.'0.1''1####Residualstandarderror:0.1299on93degreesoffreedom##MultipleR-squared:0.7063,AdjustedR-squared:0.7032##F-statistic:223.7on1and93DF,p-value:<2.2e-16所有統(tǒng)計(jì)量都是顯著的，修正R方的值是0.7032。二氧化碳水平與溫度變化高度相關(guān)。但是真的是這樣嗎？看看殘差時(shí)間序列圖：＞plot.ts(fit.lm$residuals)Time殘差時(shí)間序列圖能夠看出殘差具有明顯的模式。再看看ACF圖:acf(fit.lm$residuals)＞acf(fit.lm$residuals)：＜：＜：＜：*＜：%:v：vW77.u-n-?-j=7.U-'亠』殘差A(yù)CF圖可以看出，直到10階延遲序列，都具有顯著的自相關(guān)。對(duì)自相關(guān)性使用Durbin-Watson檢驗(yàn)：>lmtest::dwtest(fit」m)####Durbin-Watsontest####data:fit」m##DW=0.77425,p-value=4.468e-12##alternativehypothesis:trueautocorrelationisgreaterthan0這個(gè)檢驗(yàn)的原假設(shè)是“不存在自相關(guān)”，從P值可以拒絕接受原假設(shè)，即存在自相關(guān)性。2.2.2向量自回歸處理自相關(guān)的簡(jiǎn)單方式是，在相關(guān)的時(shí)間序列中加入延遲變量，使所有數(shù)據(jù)具有穩(wěn)定性。下面進(jìn)行這樣的處理，使用向量自回歸找出適當(dāng)?shù)难舆t結(jié)構(gòu)，加入我們的因果關(guān)系模型。查看數(shù)據(jù)是否平穩(wěn)，使用forecast包提供的ndiffs()函數(shù)。從它的輸出中可以找出使數(shù)據(jù)平穩(wěn)所需的最小差分次數(shù)。在這個(gè)函數(shù)中，你可以選擇3種可用的檢驗(yàn)方法，分別是：Kwiatkowski,Philips,Schmidt&Shin(KPSS))、AugmentedDickeyFuller(ADF)和Philips-Peron(PP))。使用ADF方法，它的原假設(shè)認(rèn)為數(shù)據(jù)是不平穩(wěn)的：>ndiffs(climate[,1],test="adf")##[1]1>ndiffs(climate[,2],test="adf")##[1]1在兩個(gè)序列中，一階差分即可使數(shù)據(jù)穩(wěn)定。>p_load(vars,aod)>climate.diff<-diff(climate)%>%window(start=1946)>head(climate.diff)##TimeSeries:##Start=1946##End=1951##Frequency=1##CO2Temp##194678-0.099##19471540.034##1948770.001##1949-50-0.035##1950211-0.100##19511370.121:;Ji炳詵；詵；J隸：：S：SSSS:輟W；詵;J隸€W；詵澈觀:輟鶴S；Ji:報(bào)；使用向量自回歸在信息準(zhǔn)則的基礎(chǔ)上確定最優(yōu)的延遲結(jié)構(gòu):#lag.max=12指定模型中延遲的最大數(shù)量為12lag.select<-VARseIect(cIimate.diff,lag.max=12)lag.select$selection##AIC(n)HQ(n)SC(n)FPE(n)##5115:wwwa能夠看到，AIC和FPE都選擇了5階延遲作為VAR模型的最優(yōu)結(jié)構(gòu)，HQ和SC則選擇了一階延遲。我們使用var()函數(shù)建立這個(gè)5階延遲模型:fit.5<-VAR(cIimate.diff,p=5)summary(fit.5)####VAREstimationResults:##=========================##Endogenousvariables:CO2,Temp##Deterministicvariables:const##Samplesize:63##LogLikelihood:-324.566##Rootsofthecharacteristicpolynomial:##0.84030.82150.82150.80840.80840.78430.78430.76770.62630.6263##Call:##VAR(y=climate.diff,p=5)######EstimationresultsforequationCO2:##====================================####CO2=CO2.I1+Temp.l1+CO2.I2+Temp.l2+CO2.I3+Temp.l3+CO2.I4+Temp.l4+CO2.I5+Temp.l5+const####EstimateStd.ErrortvaluePr(>|t|)##CO2.I10.292870.138242.1190.03892*##Temp.I1-61.55494150.85121-0.4080.68491##CO2.I2-0.045410.13870-0.3270.74470##Temp.I2-205.75318166.79706-1.2340.22292##CO2.I30.162320.139241.1660.24904##Temp.I3-145.88236173.15840-0.8420.40338##CO2.I4-0.091100.15722-0.5790.56478##Temp.I4-283.42320168.28093-1.6840.09813.##CO2.I50.019770.160480.1230.90243##Temp.I5-163.78328152.94657-1.0710.28918##const92.0606532.204072.8590.00611**##---##Signif.codes:0'***'0.001'**'0.01'*'0.05'.'0.1''1######Residualstandarderror:122.3on52degreesoffreedom##MultipleR-Squared:0.1741,AdjustedR-squared:0.015281.096on10and52df##F-statistic:1.096on10and52DF,p-value:0.38251.096on10and52df####F####Temp.l1-0.49247590.1240984##CO2.I20.00022980.0001141##Temp.l2-0.50638470.1372163##CO2.I3-0.00003900.0001146##Temp.l3-0.41207900.1424495##CO2.I40.0001222##Temp.l40.0423256##CO2.I5-0.0003945##Temp.l5-0.20352280.00012930.13843700.00013200.1258221-1.0200.312335p-va##EstimationresultsforequationTemp:##=============================##Temp=CO2.I1+Temp」####F####Temp.l1-0.49247590.1240984##CO2.I20.00022980.0001141##Temp.l2-0.50638470.1372163##CO2.I3-0.00003900.0001146##Temp.l3-0.41207900.1424495##CO2.I40.0001222##Temp.l40.0423256##CO2.I5-0.0003945##Temp.l5-0.20352280.00012930.13843700.00013200.1258221-1.0200.312335p-va##EstimationresultsforequationTemp:##=============================##Temp=CO2.I1+Temp」1+CO2.I2+Temp.l2+CO2.I3+Temp.l3+CO2.I4+Temp.l4+CO2.I5+Temp.l5+const##2.0140.049184*-3.6900.000537***-3.9680.000223***##EstimateStd.ErrortvaluePr(>|t|)##CO2.l1-0.00011600.0001137-0.3400.734860-2.8930.005565**0.9450.3492640.3060.761024-2.9880.004280**-1.6180.1118121.7820.080594.##const0.04720940.02649281.7820.080594.##---##Signif.codes:0'***'0.001'**'0.01'*'0.05'.'0.1''1######Residualstandarderror:0.1006on52degreesoffreedom##MultipleR-Squared:0.4509,AdjustedR-squared:0.3453##F-statistic:4.27on10and52DF,p-value:0.0002326########Covariancematrixofresiduals:##CO2Temp##CO214955.464-1.09780##Temp-1.0980.01012####Correlationmatrixofresiduals:##CO2Temp####CO21.00000-0.08923##Temp-0.089231.00000模型是顯著的(p-value:0.0002326)，修正R方為0.35。檢查序列相關(guān)性:>serial.test(fit.5,type="PT.asymptotic")####PortmanteauTest(asymptotic)####data:ResidualsofVARobjectfit.5##Chi-squared=35.912,df=44,p-value=0.8021Portmanteau檢驗(yàn)的原假設(shè)認(rèn)為自相關(guān)為0,備擇假設(shè)認(rèn)為自相關(guān)不為0°p-value=0.8021不能拒絕原假設(shè)，可以認(rèn)為殘差中不存在自相關(guān)。下面進(jìn)行格蘭杰因果關(guān)系檢驗(yàn):#建立兩個(gè)對(duì)象，一個(gè)表示x引發(fā)y,—個(gè)表示y引發(fā)xx2y<-causality(fit.5,cause="CO2")y2x<-causality(fit.5,cause="Temp")>>x2y$Granger

####GrangercausalityH0:CO2donotGranger-causeTemp####data:VARobjectfit.5##F-Test=2.2069,df1=5,df2=104,p-value=0.05908衛(wèi)羽密密密密腔密密密密旳密密密密羽密密密轉(zhuǎn)羽密密整密腔密？>y2x$Granger####GrangercausalityHO:TempdonotGranger-causeCO2####data:VARo