R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).

setwd("C:/Users/Shreya/Dropbox/Shreya/My Folder/University/Sem_8_Spring'18/Time Series/Homework 3")

library(data.table)

dt=fread('Retail_Sales_Data.csv')

plot(dt)

setnames(dt,c("Year","Month","Sales ($, Milions)"),c("year","month","sales"))
colnames(dt)        

[1] "year"  "month" "sales"

dim(dt)

[1] 194   4

Warning message:
In library(package, lib.loc = lib.loc, character.only = TRUE, logical.return = TRUE,  :
  there is no package called ‘rstream’

dt$yr_month<-seq((194))

We fit a linear tren by regressing sales against month.

lm=glm(sales~yr_month,data=dt)
summary(lm)

Call:
glm(formula = sales ~ yr_month, data = dt)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-43532   -9425      45    6457   62248  

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 140526.68    2756.66   50.98   <2e-16 ***
yr_month      1020.99      24.52   41.64   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 365717341)

    Null deviance: 7.0446e+11  on 193  degrees of freedom
Residual deviance: 7.0218e+10  on 192  degrees of freedom
AIC: 4379.7

Number of Fisher Scoring iterations: 2

Now we find the residuals of the linear fit,

res = lm$residuals

Now we will fit ARMA(m,n) usinf Fisher’s test technique to the de-trended residuals. For this we will construct a loop that does following: - Executes the Fisher’s tests in loops until it is insignificant after fitting the resquired ARMA(m,n) model. - Decided which ARMA(m,n) model to use in each loop.

plot(res)

plot(lm)

arima(res,order = c(6,0,5))

Call:
arima(x = res, order = c(6, 0, 5))

Coefficients:
         ar1      ar2     ar3      ar4     ar5      ar6      ma1     ma2
      0.7244  -0.4394  0.6445  -0.0674  0.0672  -0.1957  -0.9072  0.3886
s.e.  0.1928   0.2408  0.2265   0.2270  0.1392   0.0805   0.1970  0.2678
          ma3     ma4     ma5  intercept
      -0.8306  0.2242  0.5258   532.3564
s.e.   0.2239  0.2851  0.1984  1583.6455

sigma^2 estimated as 221215379:  log likelihood = -2145.17,  aic = 4316.35

ar0<-arima(res,order = c(1,0,0))

RSS0=sum(ar0$residuals^2)

adeq=FALSE
n=0
while(!adeq){
  n<-n+1
  print(n)
  ar1<-0
  ar1<-try(arima(res,order = c((2*n),0,((2*n)-1))))
  print(ar1)
  if(!is.integer(ar1)){
    RSS1<-sum(ar1$residuals^2)
    f_stat = ((RSS0-RSS1)/4)/(RSS1/(194- (4*n-1)))#num_data_pts - params
    if(f_stat > qf(0.95,4,(194- (4*n-1)))){
      RSS0<-RSS1
    }
    else if(f_stat > qf(0.95,4,(194- (4*n-1)))){
      print("Adequate Model found at 2n")
      print(str(2*n))
      adeq=TRUE
    }
    else{
      print("Error in Model")
      }
  }
}

[1] 1

Call:
arima(x = res, order = c((2 * n), 0, ((2 * n) - 1)))

Coefficients:
         ar1      ar2      ma1  intercept
      0.0803  -0.2389  -0.1384    75.0635
s.e.  0.2155   0.0717   0.2156   986.3519

sigma^2 estimated as 339440203:  log likelihood = -2180.69,  aic = 4371.38
[1] "Error in Model"
[1] 2

Call:
arima(x = res, order = c((2 * n), 0, ((2 * n) - 1)))

Coefficients:
         ar1      ar2      ar3      ar4      ma1     ma2     ma3  intercept
      0.2472  -0.4249  -0.5423  -0.1796  -0.4843  0.2897  0.6346   116.6167
s.e.  0.1779   0.1911   0.1772   0.1086   0.1694  0.2175  0.1782   866.8843

sigma^2 estimated as 2.53e+08:  log likelihood = -2155.87,  aic = 4329.75
[1] "Error in Model"
[1] 3

Call:
arima(x = res, order = c((2 * n), 0, ((2 * n) - 1)))

Coefficients:
         ar1      ar2     ar3      ar4     ar5      ar6      ma1     ma2      ma3
      0.7244  -0.4394  0.6445  -0.0674  0.0672  -0.1957  -0.9072  0.3886  -0.8306
s.e.  0.1928   0.2408  0.2265   0.2270  0.1392   0.0805   0.1970  0.2678   0.2239
         ma4     ma5  intercept
      0.2242  0.5258   532.3564
s.e.  0.2851  0.1984  1583.6455

sigma^2 estimated as 221215379:  log likelihood = -2145.17,  aic = 4316.35
[1] "Error in Model"
[1] 4

possible convergence problem: optim gave code = 1

Call:
arima(x = res, order = c((2 * n), 0, ((2 * n) - 1)))

Coefficients:
          ar1      ar2      ar3     ar4     ar5     ar6     ar7     ar8     ma1
      -0.2098  -0.4075  -0.5089  0.1744  0.4932  0.0920  0.7793  0.5110  0.0282
s.e.   0.0720   0.0197   0.0359  0.0150  0.0162  0.0376  0.0179  0.0732  0.0802
         ma2     ma3      ma4      ma5     ma6      ma7  intercept
      0.4896  0.9497  -0.2930  -0.0665  0.3792  -0.7648   2470.795
s.e.  0.0674  0.0911   0.1303   0.0926  0.0773   0.0796  10188.881

sigma^2 estimated as 100861266:  log likelihood = -2079.55,  aic = 4193.1
[1] "Error in Model"
[1] 5

possible convergence problem: optim gave code = 1

Call:
arima(x = res, order = c((2 * n), 0, ((2 * n) - 1)))

Coefficients:
         ar1      ar2     ar3      ar4     ar5     ar6     ar7      ar8     ar9
      0.1088  -0.6904  0.5628  -0.1811  0.9346  0.0282  0.8238  -0.2409  0.3439
s.e.  0.0504   0.0428  0.0627   0.0443  0.0426  0.0477  0.0490   0.0646  0.0459
         ar10      ma1     ma2      ma3     ma4      ma5     ma6      ma7      ma8
      -0.7578  -0.0891  1.3061  -0.4691  0.9272  -1.0046  0.3510  -1.2622  -0.0997
s.e.   0.0475   0.0817  0.0671   0.1064  0.1023   0.0934  0.0953   0.0704   0.0691
          ma9  intercept
      -0.6577     0.7247
s.e.   0.0710   836.6481

sigma^2 estimated as 46764684:  log likelihood = -2003.74,  aic = 4049.48
[1] "Error in Model"
[1] 6

possible convergence problem: optim gave code = 1

Call:
arima(x = res, order = c((2 * n), 0, ((2 * n) - 1)))

Coefficients:
          ar1      ar2      ar3     ar4      ar5      ar6      ar7     ar8      ar9
      -0.0355  -0.0180  -0.0467  -0.051  -0.0093  -0.0368  -0.0305  0.0185  -0.0339
s.e.   0.0195   0.0202   0.0211   0.020   0.0208   0.0208   0.0208  0.0207   0.0211
         ar10    ar11    ar12     ma1     ma2     ma3     ma4     ma5     ma6
      -0.0170  0.0076  0.9535  0.1774  0.1056  0.4706  0.2108  0.3350  0.4154
s.e.   0.0205  0.0199  0.0190  0.0787  0.0813  0.0707  0.0740  0.0727  0.0708
         ma7     ma8     ma9    ma10     ma11  intercept
      0.3457  0.1743  0.6116  0.2118  -0.0126   1992.580
s.e.  0.0753  0.0850  0.0612  0.0872   0.0900   3509.994

sigma^2 estimated as 17757898:  log likelihood = -1916.23,  aic = 3882.46
[1] "Error in Model"
[1] 7

possible convergence problem: optim gave code = 1

Call:
arima(x = res, order = c((2 * n), 0, ((2 * n) - 1)))

Coefficients:

NaNs produced

         ar1     ar2      ar3      ar4     ar5      ar6      ar7     ar8      ar9
      0.2559  0.6383  -0.0207  -0.0328  0.0245  -0.0012  -0.0241  0.0437  -0.0216
s.e.     NaN     NaN   0.0089   0.0190  0.0187   0.0246   0.0232  0.0056   0.0282
         ar10    ar11    ar12     ar13     ar14      ma1      ma2     ma3     ma4
      -0.0299  0.0298  0.9586  -0.2907  -0.6242  -0.1486  -0.6158  0.3598  0.0255
s.e.   0.0169     NaN  0.0175      NaN      NaN      NaN      NaN  0.0868     NaN
          ma5     ma6     ma7      ma8     ma9     ma10     ma11     ma12    ma13
      -0.0081  0.1895  0.0448  -0.1934  0.3509  -0.0756  -0.4866  -0.0568  0.0652
s.e.      NaN     NaN  0.0938   0.0441  0.0121      NaN      NaN      NaN  0.0895
      intercept
        580.693
s.e.   1560.413

sigma^2 estimated as 17234028:  log likelihood = -1914.44,  aic = 3886.87
[1] "Error in Model"
[1] 8

Error in optim(init[mask], armafn, method = optim.method, hessian = TRUE,  : 
  non-finite finite-difference value [16]

[1] "Error in optim(init[mask], armafn, method = optim.method, hessian = TRUE,  : \n  non-finite finite-difference value [16]\n"
attr(,"class")
[1] "try-error"
attr(,"condition")
<simpleError in optim(init[mask], armafn, method = optim.method, hessian = TRUE,     control = optim.control, trans = as.logical(transform.pars)): non-finite finite-difference value [16]>

Error: $ operator is invalid for atomic vectors