## Regression with time series error term. #Look at the example in the book. TInt1year=read.csv("1yearTreasury1962to2009.csv",header=T) ### These are weekly data given in percents TInt1year[1,] TInt1year=read.csv("1yearTreasury1962to2009.csv",header=T,skip=7) TInt1year=read.csv("1yearTreasury1962to2009.csv",header=T,skip=7)[,2] TInt3year=read.csv("3yearTreasury1962to2009.csv",header=T,skip=7) [,2] plot(TInt1year,type="l",col="green") lines(TInt3year,col="red") ## If you want some semblance of the years on the x axis plot(seq(1962,2009.62,by=1/52),TInt1year,type="l",col="green") lines(seq(1962,2009.62,by=1/52),TInt3year,col="red") ## Of course you could use the data stuff date.Tbonds=strptime(read.csv("1yearTreasury1962to2009.csv",header=T,skip=7)[,1],"%m/%d/%Y") plot(date.Tbonds,TInt1year,type="l",col="green") lines(date.Tbonds,TInt3year,col="red") plot(TInt1year,TInt3year,lty=15,pch=20) ## a linear relation seems very appropriate plot(diff(TInt1year),diff(TInt3year),lty=15,pch=20) ## changes in the rates (Figure 2.18 in the book) #So let us try regression: model1=lm(TInt3year~TInt1year) summary(model1) ## Seemingly very good model with very high R squared. plot(model1) ## regular diagnostic plots - do not exhibit much except for the normal plot ## Now let us plot residuals in order and test for white noise plot(model1$residuals,type="l") acf(model1$residuals) ## clearly not white noise Box.test(model1$residuals,lag=10,type="Ljung") ## recall: null hypothesis= data is not correlated for the first x lags ## unit root testing: library(fUnitRoots) ar(TInt1year) adfTest(TInt1year,lag=30) ## seems there is a unit root #So consider differences instead c1=diff(TInt1year) c3=diff(TInt3year) model2= lm(c3~c1) summary(model2) Call: lm(formula = c3 ~ c1) Residuals: Min 1Q Median 3Q Max -0.424922 -0.036154 -0.001413 0.034495 0.489375 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.0001142 0.0013933 0.082 0.935 c1 0.7922140 0.0073760 107.404 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.06933 on 2474 degrees of freedom Multiple R-squared: 0.8234, Adjusted R-squared: 0.8233 F-statistic: 1.154e+04 on 1 and 2474 DF, p-value: < 2.2e-16 model2= lm(c3~c1-1) ## the model without the intercept summary(model2) plot(model2) ## actually these plots look worse than the previous ones plot(model2$residuals,type="l") acf(model2$residuals) ## this plot now looks like a tipical time series - clearly not white noise though which is why we get those types of plots before ## Now we can use the methods before for the residual series and build a time series model for them. library(tseries) eacf.tables(model2$residuals) [1] "EACF table" [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [1,] 0.17476776 -0.001469869 -0.002503988 0.06406689 -3.079153e-02 -0.0934975974 -0.07700085 0.0019587977 0.0007077368 -0.053309392 -0.0184702810 0.0417679567 0.0671271769 -0.01119983 0.023948796 [2,] 0.18739204 -0.105122902 0.001054763 0.04517043 -2.574459e-02 -0.0482399456 -0.07840888 0.0354282855 -0.0003326228 -0.053470463 -0.0201770655 0.0242407047 0.0682596238 0.02126555 0.022357684 [3,] 0.10626307 0.010020378 0.005477349 0.01417748 -4.880270e-02 0.0039617587 -0.07806416 0.0001547369 -0.0055617805 -0.028239692 0.0112291803 -0.0041718476 0.0673959485 -0.04661191 0.024359512 [4,] -0.05093700 0.357248106 -0.009230458 0.03614059 -4.841580e-02 0.0001042357 -0.04510199 -0.0113304962 -0.0070213087 -0.044145873 -0.0043081928 0.0090221368 0.0769128544 -0.04376635 0.027422642 [5,] 0.49609041 0.371103020 0.024379910 0.22669888 -8.080847e-02 -0.0382817958 -0.04619610 -0.0012547288 -0.0051840746 -0.043102319 -0.0036066200 0.0161510636 0.0730534597 0.02917100 0.014379088 [6,] -0.47064703 0.204120025 0.129791074 -0.26058314 -9.184419e-03 0.0337549741 -0.04770678 0.0057997395 0.0023224437 -0.004217892 0.0004798517 -0.0154956991 0.0695500913 -0.01061674 0.021969693 [7,] -0.44673379 0.260631220 0.193792865 -0.25165623 -1.540023e-02 0.2715201280 0.05474164 -0.0133050184 0.0122146265 0.008963857 -0.0109948991 -0.0057883807 0.0021324221 -0.03726367 0.026878617 [8,] 0.30525230 0.064890458 -0.273124670 0.02771996 -4.709124e-02 0.3771704909 0.12387807 0.0177056577 0.0120704189 0.006395117 -0.0182178490 -0.0168927764 -0.0047652512 -0.04948188 0.027603809 [9,] 0.02822761 0.035305353 -0.273983335 -0.21041304 -4.332714e-05 0.3018449123 0.17259360 -0.0795544478 0.0243662200 0.001164909 -0.0241431445 -0.0040502488 -0.0006889595 -0.05043751 -0.019087516 [10,] -0.00554319 0.315720499 -0.349265176 -0.04410145 2.334735e-02 0.1742203334 0.03000833 0.1492582121 0.0252440751 0.001596173 -0.0055238359 -0.0069978439 -0.0142929866 -0.03876004 0.032816813 [11,] 0.03562720 0.465462226 0.372025379 0.07415333 1.837405e-01 0.1524898760 -0.01189708 0.1758226603 -0.0687158441 0.031752449 0.0058311557 -0.0003655077 -0.0111025732 -0.03664001 0.031972124 [12,] -0.04689873 0.177330271 0.405435685 -0.09737710 1.458915e-01 0.1718069500 -0.02763217 0.0816165405 -0.1396956268 0.084271576 -0.0042185091 0.0056913339 -0.0112002302 -0.02951228 0.030653927 [13,] -0.46767056 0.231061275 0.232343816 -0.26717093 2.355394e-01 0.2396091810 0.13524463 -0.0053118288 -0.1430050951 0.038415702 0.0499742460 0.0500826947 -0.0017026763 -0.04487658 0.015853952 [14,] 0.43959865 -0.160157478 -0.287098214 -0.08942529 5.123660e-02 0.2796924339 0.06355757 -0.0164947052 -0.1790050723 0.149276810 -0.1007372255 0.0816020063 0.0354260657 -0.04316881 -0.006352550 [15,] 0.48964371 -0.338530828 -0.280733011 -0.14858955 7.128456e-03 0.1382256589 0.25374250 0.2296207030 -0.1056252625 -0.132454457 -0.1466799244 -0.1480782718 -0.2404290136 -0.05434749 -0.002152781 [1] " " [1] "Simplified EACF: 2 denotes significance" [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [1,] 2 0 0 2 0 2 2 0 0 2 0 2 2 0 0 [2,] 2 2 0 2 0 2 2 0 0 2 0 0 2 0 0 [3,] 2 0 0 0 2 0 2 0 0 0 0 0 2 2 0 [4,] 2 2 0 0 2 0 2 0 0 2 0 0 2 2 0 [5,] 2 2 0 2 2 0 2 0 0 2 0 0 2 0 0 [6,] 2 2 2 2 0 0 2 0 0 0 0 0 2 0 0 [7,] 2 2 2 2 0 2 2 0 0 0 0 0 0 0 0 [8,] 2 2 2 0 2 2 2 0 0 0 0 0 0 2 0 [9,] 0 0 2 2 0 2 2 2 0 0 0 0 0 2 0 [10,] 0 2 2 2 0 2 0 2 0 0 0 0 0 0 0 [11,] 0 2 2 2 2 2 0 2 2 0 0 0 0 0 0 [12,] 2 2 2 2 2 2 0 2 2 2 0 0 0 0 0 [13,] 2 2 2 2 2 2 2 0 2 0 2 2 0 2 0 [14,] 2 2 2 2 2 2 2 0 2 2 2 2 0 2 0 [15,] 2 2 2 2 0 2 2 2 2 2 2 2 2 2 0 Warning message: In optim(coef, err, gr = NULL, hessian = TRUE, ...) : one-diml optimization by Nelder-Mead is unreliable: use optimize ##Note I'll also looka at the other included eacf function: library(TSA) eacf(model2$residuals) AR/MA 0 1 2 3 4 5 6 7 8 9 10 11 12 13 0 x o o x o x x o o x o x x o 1 x x o x o x x o o x o o x o 2 x o o o x o x o o o o o x x 3 x x o o x o x o o x o o x x 4 x x o x x o x o o x o o x o 5 x x x x o o x o o o o o x o 6 x x x x o x x o o o o o o o 7 x x x o x x x o o o o o o x ## Perhaps MA(1) is appropriate modelresid=arima(model2$residuals,order=c(0,0,1),include.mean=FALSE) modelresid Call: arima(x = model2$residuals, order = c(0, 0, 1), include.mean = FALSE) Coefficients: ma1 0.1796 s.e. 0.0196 sigma^2 estimated as 0.004651: log likelihood = 3135.6, aic = -6269.19 arima(model2$residuals,order=c(0,0,1),include.mean=T) ## In fact the mean is zero as this shows ## So the model is: ## c3= 0.7922140 c1 + e where e_t=a_t+ 0.1796 a_{t-1} and a_t a white noise ## rewriting the model in terms of r1 and r3 and using the fact that the intercept is zero we get: ## r3_t=r3_{t-1}+0.792*(r1_t-r1_{t-1}) + a_t+0.1796 a_{t-1}