This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

summary

```

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot. #import the data on the growth rate of GDP, convert it into time series ts object

#download data

library("Quandl")
## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
#(1) data identification 
 
rcpg <- Quandl("FRED/PCECC96", type="ts")

summary(rcpg) 
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1199    2212    4050    5004    7465   11320
plot(rcpg)

str(rcpg)
##  Time-Series [1:276] from 1947 to 2016: 1199 1219 1223 1224 1230 ...
head(rcpg)
## [1] 1199.4 1219.3 1223.3 1223.6 1229.8 1244.1
tail(rcpg)
## [1] 10918.6 11033.3 11081.2 11178.9 11262.4 11322.5
# transformed data

g=log(rcpg) 
 
 Y=diff(g)
 Y
##               Qtr1          Qtr2          Qtr3          Qtr4
## 1947                0.0164554918  0.0032752015  0.0002452082
## 1948  0.0050542212  0.0115608224  0.0014457834  0.0079146594
## 1949  0.0016708442  0.0151482187  0.0022681954  0.0145806833
## 1950  0.0164939598  0.0163008304  0.0500631214 -0.0306564987
## 1951  0.0239009175 -0.0286623080  0.0115346585  0.0058621481
## 1952  0.0022344767  0.0193945473  0.0048607049  0.0347374636
## 1953  0.0118081059  0.0060188773 -0.0023363717 -0.0066383051
## 1954  0.0036263555  0.0129867631  0.0134733069  0.0208031276
## 1955  0.0221943289  0.0189549778  0.0122581196  0.0124689144
## 1956  0.0016149778  0.0033412919  0.0022609637  0.0136939928
## 1957  0.0068355362  0.0017452011  0.0078736006  0.0004612280
## 1958 -0.0138133363  0.0081485813  0.0163287708  0.0133690316
## 1959  0.0180124764  0.0153018040  0.0103958920  0.0011305215
## 1960  0.0094782364  0.0126047639 -0.0040078100  0.0012673604
## 1961 -0.0002639010  0.0147762709  0.0048772991  0.0198330738
## 1962  0.0105997358  0.0121764067  0.0081016103  0.0140702893
## 1963  0.0067206829  0.0095437504  0.0135585655  0.0082999703
## 1964  0.0195150702  0.0174328003  0.0181952933  0.0028243620
## 1965  0.0220536286  0.0110176237  0.0170778129  0.0276945864
## 1966  0.0146502013  0.0025278376  0.0113166652  0.0041123031
## 1967  0.0057838882  0.0135605866  0.0050959459  0.0060657434
## 1968  0.0234904418  0.0151383122  0.0186062060  0.0045080048
## 1969  0.0111285209  0.0063848117  0.0048323757  0.0079558818
## 1970  0.0061598622  0.0046037569  0.0087679142 -0.0027424514
## 1971  0.0189737076  0.0091199291  0.0079453885  0.0164858747
## 1972  0.0131372218  0.0189147495  0.0153071400  0.0231829471
## 1973  0.0181073916 -0.0004173996  0.0035423556 -0.0029163216
## 1974 -0.0087702794  0.0035113555  0.0040959770 -0.0147580863
## 1975  0.0083225762  0.0165583461  0.0141942029  0.0105437932
## 1976  0.0197998024  0.0091391607  0.0105532326  0.0129889825
## 1977  0.0113637587  0.0054994073  0.0094985262  0.0149599724
## 1978  0.0057549599  0.0211120960  0.0041796279  0.0079792710
## 1979  0.0050584598 -0.0005775339  0.0097730010  0.0026826982
## 1980 -0.0015391875 -0.0227411019  0.0107188123  0.0131644941
## 1981  0.0052472770 -0.0001728203  0.0040165150 -0.0075287620
## 1982  0.0065938376  0.0036854173  0.0076954464  0.0180876006
## 1983  0.0096802954  0.0195946831  0.0173949442  0.0156389332
## 1984  0.0084526442  0.0141504495  0.0076546608  0.0131380062
## 1985  0.0168655320  0.0093436989  0.0190256675  0.0025656565
## 1986  0.0084304279  0.0111177390  0.0179499406  0.0063768446
## 1987  0.0001569397  0.0137731686  0.0115225712  0.0021016439
## 1988  0.0176316026  0.0073053714  0.0085083233  0.0113789838
## 1989  0.0046064152  0.0046937808  0.0098950145  0.0043942767
## 1990  0.0085543417  0.0031230451  0.0040261313 -0.0075910716
## 1991 -0.0034535008  0.0083564223  0.0048439843 -0.0002626579
## 1992  0.0185996999  0.0068354371  0.0107661214  0.0118372396
## 1993  0.0037817936  0.0089392729  0.0109813766  0.0088122025
## 1994  0.0114064791  0.0077176225  0.0077214364  0.0107014805
## 1995  0.0026420505  0.0089311141  0.0091264702  0.0070025425
## 1996  0.0092360967  0.0107997887  0.0060055799  0.0078298122
## 1997  0.0104949253  0.0045219855  0.0169654264  0.0118062797
## 1998  0.0102693626  0.0175069399  0.0130596977  0.0145888615
## 1999  0.0094821191  0.0146971415  0.0112921436  0.0145748895
## 2000  0.0151106759  0.0095508878  0.0096797648  0.0088629738
## 2001  0.0042159086  0.0025689982  0.0036381084  0.0151630321
## 2002  0.0029958257  0.0050899032  0.0069667241  0.0053634306
## 2003  0.0043826169  0.0110719086  0.0146377882  0.0077334046
## 2004  0.0096768535  0.0064760607  0.0095117167  0.0102041702
## 2005  0.0076172556  0.0108136588  0.0077186494  0.0037591485
## 2006  0.0111522822  0.0053100554  0.0058208747  0.0101434389
## 2007  0.0052486184  0.0033874119  0.0044391875  0.0012505584
## 2008 -0.0020652548  0.0016783443 -0.0072499446 -0.0121068558
## 2009 -0.0034354370 -0.0045174364  0.0060491332 -0.0001115014
## 2010  0.0053481740  0.0081040406  0.0064501881  0.0101833874
## 2011  0.0050041315  0.0020141978  0.0043372599  0.0033593895
## 2012  0.0060108995  0.0016942956  0.0026416034  0.0027877186
## 2013  0.0061989384  0.0034262698  0.0042168623  0.0085959203
## 2014  0.0032033499  0.0094288667  0.0084893979  0.0104502156
## 2015  0.0043320061  0.0087780942  0.0074416710  0.0053221526
 plot(Y)

# examine ACF and PACF of the transformed data to determine plausible models to be estimated 
#acf(Y),we find the order of AR by two ways 1- PACF here we find the order is 2 from PACF plot, #2-information cariteria AICtransformed here we prefer AR(4) based on AIC.we also find data not #correlated at AR(1) but at AR(2) and AR(3) correlated data are correlated based on L-b test
#from acf we choose ma(1)
#from pacf we choose ar(2)
#try Ar(2),Ar(3),Ar(4)
#we use AIC to choose order p of Ar(p) we find the order is 4 Ar(4)
#use BIC to choose p of Ar(p)
#AIC and BIC prefer Ar(4)

acf(as.data.frame(Y),type='correlation',lag=24)

 pacf(Y)

 A=arima(Y,order=c(1,0,0),method="ML")
 A
## 
## Call:
## arima(x = Y, order = c(1, 0, 0), method = "ML")
## 
## Coefficients:
##          ar1  intercept
##       0.0893     0.0082
## s.e.  0.0601     0.0005
## 
## sigma^2 estimated as 6.649e-05:  log likelihood = 932.32,  aic = -1858.64
 str(A)
## List of 14
##  $ coef     : Named num [1:2] 0.08933 0.00816
##   ..- attr(*, "names")= chr [1:2] "ar1" "intercept"
##  $ sigma2   : num 6.65e-05
##  $ var.coef : num [1:2, 1:2] 3.61e-03 7.69e-08 7.69e-08 2.95e-07
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2] "ar1" "intercept"
##   .. ..$ : chr [1:2] "ar1" "intercept"
##  $ mask     : logi [1:2] TRUE TRUE
##  $ loglik   : num 932
##  $ aic      : num -1859
##  $ arma     : int [1:7] 1 0 0 0 4 0 0
##  $ residuals: Time-Series [1:275] from 1947 to 2016: 0.00826 -0.00563 -0.00748 -0.0024 0.00367 ...
##  $ call     : language arima(x = Y, order = c(1, 0, 0), method = "ML")
##  $ series   : chr "Y"
##  $ code     : int 0
##  $ n.cond   : num 0
##  $ nobs     : int 275
##  $ model    :List of 10
##   ..$ phi  : num 0.0893
##   ..$ theta: num(0) 
##   ..$ Delta: num(0) 
##   ..$ Z    : num 1
##   ..$ a    : num -0.00284
##   ..$ P    : num [1, 1] 0
##   ..$ T    : num [1, 1] 0.0893
##   ..$ V    : num [1, 1] 1
##   ..$ h    : num 0
##   ..$ Pn   : num [1, 1] 1
##  - attr(*, "class")= chr "Arima"
 B=arima(Y,order=c(0,0,2),method="ML")
 B
## 
## Call:
## arima(x = Y, order = c(0, 0, 2), method = "ML")
## 
## Coefficients:
##          ma1     ma2  intercept
##       0.0267  0.3660     0.0082
## s.e.  0.0567  0.0586     0.0006
## 
## sigma^2 estimated as 5.889e-05:  log likelihood = 948.88,  aic = -1889.77
 str(B)
## List of 14
##  $ coef     : Named num [1:3] 0.02674 0.36602 0.00817
##   ..- attr(*, "names")= chr [1:3] "ma1" "ma2" "intercept"
##  $ sigma2   : num 5.89e-05
##  $ var.coef : num [1:3, 1:3] 3.22e-03 -5.08e-04 1.28e-07 -5.08e-04 3.43e-03 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:3] "ma1" "ma2" "intercept"
##   .. ..$ : chr [1:3] "ma1" "ma2" "intercept"
##  $ mask     : logi [1:3] TRUE TRUE TRUE
##  $ loglik   : num 949
##  $ aic      : num -1890
##  $ arma     : int [1:7] 0 2 0 0 4 0 0
##  $ residuals: Time-Series [1:275] from 1947 to 2016: 0.00778 -0.00485 -0.0104 -0.00115 0.00719 ...
##  $ call     : language arima(x = Y, order = c(0, 0, 2), method = "ML")
##  $ series   : chr "Y"
##  $ code     : int 0
##  $ n.cond   : num 0
##  $ nobs     : int 275
##  $ model    :List of 10
##   ..$ phi  : num(0) 
##   ..$ theta: num [1:2] 0.0267 0.366
##   ..$ Delta: num(0) 
##   ..$ Z    : num [1:3] 1 0 0
##   ..$ a    : num [1:3] -0.002848 0.000258 -0.001051
##   ..$ P    : num [1:3, 1:3] 0 0 0 0 0 0 0 0 0
##   ..$ T    : num [1:3, 1:3] 0 0 0 1 0 0 0 1 0
##   ..$ V    : num [1:3, 1:3] 1 0.026737 0.36602 0.026737 0.000715 ...
##   ..$ h    : num 0
##   ..$ Pn   : num [1:3, 1:3] 1 0.026737 0.36602 0.026737 0.000715 ...
##  - attr(*, "class")= chr "Arima"

```r #Check model A for adequacy, some serial correlation not significant ACF and low p-values for Ljunb-Box test at many lags

tsdiag(A,gof.lag=24) ```

```r # try AR(2)

A2 <- arima(Y, order=c(2,0,0))

A2 ```

## ## Call: ## arima(x = Y, order = c(2, 0, 0)) ## ## Coefficients: ## ar1 ar2 intercept ## 0.0599 0.3188 0.0082 ## s.e. 0.0571 0.0570 0.0007 ## ## sigma^2 estimated as 5.968e-05: log likelihood = 947.08, aic = -1886.16

```r # Check model A for adequacy, some serial correlation not significant ACF and low p-values for Ljunb-Box test at many lags. but this better than AR(1)

tsdiag(A2,gof.lag=24) ```

```r #try AR(3)

A3 <- arima(Y, order=c(3,0,0))

A3 ```

## ## Call: ## arima(x = Y, order = c(3, 0, 0)) ## ## Coefficients: ## ar1 ar2 ar3 intercept ## 0.0545 0.3178 0.0165 0.0082 ## s.e. 0.0604 0.0571 0.0603 0.0008 ## ## sigma^2 estimated as 5.966e-05: log likelihood = 947.12, aic = -1884.23

#Check model A for adequacy,some serial correlation not significant ACF and low p-values for Ljunb-Box test at many lags. but this better than AR(1)  

 tsdiag(A3,gof.lag=24)

 A3$coef/sqrt(diag(A3$var.coef))
##        ar1        ar2        ar3  intercept 
##  0.9031198  5.5618712  0.2738954 10.7257884
 (1-pnorm(abs(A3$coef)/sqrt(diag(A3$var.coef))))*2
##          ar1          ar2          ar3    intercept 
## 3.664623e-01 2.668972e-08 7.841650e-01 0.000000e+00
#try AR(4)
 
A4 <- arima(Y, order=c(4,0,0))
A4
## 
## Call:
## arima(x = Y, order = c(4, 0, 0))
## 
## Coefficients:
##          ar1     ar2     ar3      ar4  intercept
##       0.0570  0.3647  0.0244  -0.1444     0.0082
## s.e.  0.0597  0.0598  0.0598   0.0596     0.0007
## 
## sigma^2 estimated as 5.84e-05:  log likelihood = 950.02,  aic = -1888.03
#Check model A4 for adequacy,all serial correlation  significant ACF the problems with lags in AR(1), AR(2), and AR(3)are gone.

 tsdiag(A4,gof.lag=24)

 A4$coef/sqrt(diag(A4$var.coef))
##        ar1        ar2        ar3        ar4  intercept 
##  0.9554484  6.0975802  0.4076225 -2.4215870 12.3531993
 (1-pnorm(abs(A4$coef)/sqrt(diag(A4$var.coef))))*2
##          ar1          ar2          ar3          ar4    intercept 
## 3.393510e-01 1.076861e-09 6.835508e-01 1.545290e-02 0.000000e+00
 m <- ar(Y,method="mle")
 
 str(m)
## List of 14
##  $ order       : int 4
##  $ ar          : num [1:4] 0.057 0.3647 0.0244 -0.1444
##  $ var.pred    : num 5.84e-05
##  $ x.mean      : num 0.00816
##  $ aic         : Named num [1:13] 29.59 29.39 1.87 3.81 0 ...
##   ..- attr(*, "names")= chr [1:13] "0" "1" "2" "3" ...
##  $ n.used      : int 275
##  $ order.max   : num 12
##  $ partialacf  : NULL
##  $ resid       : Time-Series [1:275] from 1947 to 2016: NA NA NA NA 0.00778 ...
##  $ method      : chr "MLE"
##  $ series      : chr "Y"
##  $ frequency   : num 4
##  $ call        : language ar(x = Y, method = "mle")
##  $ asy.var.coef: num [1:4, 1:4] 3.56e-03 -1.97e-04 -1.13e-03 -5.94e-05 -1.97e-04 ...
##  - attr(*, "class")= chr "ar"
 m
## 
## Call:
## ar(x = Y, method = "mle")
## 
## Coefficients:
##       1        2        3        4  
##  0.0570   0.3647   0.0244  -0.1444  
## 
## Order selected 4  sigma^2 estimated as  5.84e-05
 m$order
## [1] 4
 m$aic
##          0          1          2          3          4          5 
## 29.5898061 29.3905747  1.8742608  3.8053769  0.0000000  1.9999836 
##          6          7          8          9         10         11 
##  2.8643601  1.8839187  1.6397094  1.8639941  0.5695805  2.5608319 
##         12 
##  1.1372523
 BIC(A2)
## [1] -1871.691
 BIC(A3)
## [1] -1866.149
 BIC(A4)
## [1] -1866.332
 A2.LB.lag12 <- Box.test(A2$residuals, lag=12, type="Ljung")
 
 str(A2.LB.lag12)
## List of 5
##  $ statistic: Named num 19.3
##   ..- attr(*, "names")= chr "X-squared"
##  $ parameter: Named num 12
##   ..- attr(*, "names")= chr "df"
##  $ p.value  : num 0.0818
##  $ method   : chr "Box-Ljung test"
##  $ data.name: chr "A2$residuals"
##  - attr(*, "class")= chr "htest"
 A2.LB.lag12
## 
##  Box-Ljung test
## 
## data:  A2$residuals
## X-squared = 19.286, df = 12, p-value = 0.08185
 1-pchisq(A2.LB.lag12$statistic,10)
##  X-squared 
## 0.03677295
 A2.LB.lag16 <- Box.test(A2$residuals, lag=16, type="Ljung")
 
 A2.LB.lag16
## 
##  Box-Ljung test
## 
## data:  A2$residuals
## X-squared = 24.683, df = 16, p-value = 0.07561
 1-pchisq(A2.LB.lag16$statistic,14)
##  X-squared 
## 0.03782697