Perbandingan Model Koyck, Regression with Distributed Lag 1 & 2, dan Autoregressive

Dinda Khamila Nurfatimah

2023-09-12

Packages

## Loading required package: nardl

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

## Loading required package: dynlm

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## 
## Attaching package: 'MLmetrics'

## The following object is masked from 'package:dLagM':
## 
##     MAPE

## The following object is masked from 'package:base':
## 
##     Recall

## Loading required package: carData

Impor Data

library(rio)
data <- import("https://raw.githubusercontent.com/DindaKhamila/mpdw/main/Pertemuan%203/data3.csv")
str(data)

## 'data.frame':    72 obs. of  2 variables:
##  $ AQI: num  30.2 28.2 26.6 25 26 26 26 24 23 24 ...
##  $ CO : num  199 198 199 202 205 ...

data

##     AQI       CO
## 1  30.2 198.6027
## 2  28.2 197.6013
## 3  26.6 198.6027
## 4  25.0 201.9405
## 5  26.0 205.2784
## 6  26.0 208.6163
## 7  26.0 208.6163
## 8  24.0 205.2784
## 9  23.0 200.2716
## 10 24.0 196.9338
## 11 25.0 198.6027
## 12 26.0 200.2716
## 13 27.0 203.6095
## 14 29.0 205.2784
## 15 29.0 205.2784
## 16 29.0 205.2784
## 17 29.0 203.6095
## 18 28.0 201.9405
## 19 27.0 200.2716
## 20 27.0 201.9405
## 21 27.0 201.9405
## 22 27.0 201.9405
## 23 27.0 200.2716
## 24 27.0 200.2716
## 25 27.0 200.2716
## 26 28.0 201.9405
## 27 29.0 203.6095
## 28 31.0 205.2784
## 29 31.0 206.9473
## 30 32.0 206.9473
## 31 32.0 205.2784
## 32 31.0 206.9473
## 33 31.0 206.9473
## 34 30.0 208.6163
## 35 30.0 206.9473
## 36 29.0 205.2784
## 37 27.0 203.6095
## 38 26.0 200.2716
## 39 25.0 200.2716
## 40 25.0 198.6027
## 41 24.0 196.9338
## 42 24.0 195.2648
## 43 25.0 195.2648
## 44 25.0 195.2648
## 45 25.0 196.9338
## 46 26.0 198.6027
## 47 26.0 198.6027
## 48 27.0 198.6027
## 49 27.0 198.6027
## 50 27.0 198.6027
## 51 27.0 200.2716
## 52 28.0 200.2716
## 53 28.0 200.2716
## 54 27.0 200.2716
## 55 27.0 198.6027
## 56 26.0 198.6027
## 57 25.0 198.6027
## 58 23.0 198.6027
## 59 22.0 196.9338
## 60 21.0 193.5959
## 61 20.0 193.5959
## 62 19.0 191.9270
## 63 19.0 191.9270
## 64 20.0 191.9270
## 65 21.0 191.9270
## 66 21.0 191.9270
## 67 21.0 191.9270
## 68 22.0 193.5959
## 69 24.0 195.2648
## 70 26.0 196.9338
## 71 27.0 198.6027
## 72 28.0 198.6027

Pembagian Data

#SPLIT DATA
train<-data[1:57,]
test<-data[58:72,]

#data time series
train.ts<-ts(train)
test.ts<-ts(test)
data.ts<-ts(data)

Model Koyck

Model Koyck didasarkan pada asumsi bahwa semakin jauh jarak lag peubah independen dari periode sekarang maka semakin kecil pengaruh peubah lag terhadap peubah dependen.

Koyck mengusulkan suatu metode untuk menduga model dinamis distributed lag dengan mengasumsikan bahwa semua koefisien \(\beta\) mempunyai tanda sama.

Model kyock merupakan jenis paling umum dari model infinite distributed lag dan juga dikenal sebagai geometric lag

\[ y_t=a(1-\lambda)+\beta_0X_t+\beta_1Z_t+\lambda Y_{t-1}+V_t \]

dengan \[V_t=u_t-\lambda u_{t-1}\]

Pemodelan

Pemodelan model Koyck dengan R dapat menggunakan dLagM::koyckDlm() . Fungsi umum dari koyckDlm adalah sebagai berikut.

koyckDlm(x , y , intercept)

Fungsi koyckDlm() akan menerapkan model lag terdistribusi dengan transformasi Koyck satu prediktor. Nilai x dan y tidak perlu sebagai objek time series (ts). intercept dapat dibuat TRUE untuk memasukkan intersep ke dalam model.

#MODEL KOYCK
model.koyck <- koyckDlm(x = train$CO, y = train$AQI)
summary(model.koyck)

## 
## Call:
## "Y ~ (Intercept) + Y.1 + X.t"
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.0320 -0.6498  0.0663  0.5993  2.2743 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.8427209  9.6005572   0.296    0.768    
## Y.1          0.8979076  0.0828445  10.838 4.67e-15 ***
## X.t         -0.0007616  0.0549301  -0.014    0.989    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9314 on 53 degrees of freedom
## Multiple R-Squared: 0.8177,  Adjusted R-squared: 0.8108 
## Wald test: 118.9 on 2 and 53 DF,  p-value: < 2.2e-16 
## 
## Diagnostic tests:
## NULL
## 
##                            alpha          beta       phi
## Geometric coefficients:  27.8446 -0.0007616174 0.8979076

AIC(model.koyck)

## [1] 155.8773

BIC(model.koyck)

## [1] 163.9787

Dari hasil tersebut, didapat bahwa peubah \(y_{t-1} (4.67e-15)\) memiliki nilai \(P-Value<0.05\). Hal ini menunjukkan bahwa peubah \(y_{t-1}\) berpengaruh signifikan terhadap \(y\). Adapun model keseluruhannya adalah sebagai berikut

\[ \hat{Y_t}=2.8427209-0.0007616X_t+0.8979076Y_{t-1} \]

Peramalan dan Akurasi

Berikut adalah hasil peramalan y untuk 5 periode kedepan menggunakan model koyck

fore.koyck <- forecast(model = model.koyck, x=test$CO, h=15)
fore.koyck

## $forecasts
##  [1] 25.13915 25.26537 25.38124 25.48529 25.57998 25.66500 25.74135 25.80990
##  [9] 25.87145 25.92672 25.97507 26.01722 26.05379 26.08536 26.11371
## 
## $call
## forecast.koyckDlm(model = model.koyck, x = test$CO, h = 15)
## 
## attr(,"class")
## [1] "forecast.koyckDlm" "dLagM"

mape.koyck <- MAPE(fore.koyck$forecasts, test$AQI)
#akurasi data training
GoF(model.koyck)

##              n      MAE          MPE      MAPE      sMAPE     MASE       MSE
## model.koyck 56 0.692456 -0.001096694 0.0256365 0.02556864 1.081963 0.8210179
##                  MRAE    GMRAE
## model.koyck 787601328 12509.62

Pada perhitungan keakuratan model menggunakan metode Koyck didapatkan nilai MAPE 2,56%. Nilai akurasi model ini kurang dari 10% sehingga dapat dikategorikan sangat baik.

Regression with Distributed Lag

Pemodelan model Regression with Distributed Lag dengan R dapat menggunakan dLagM::dlm() . Fungsi umum dari dlm adalah sebagai berikut.

dlm(formula , data , x , y , q , remove )

Fungsi dlm() akan menerapkan model lag terdistribusi dengan satu atau lebih prediktor. Nilai x dan y tidak perlu sebagai objek time series (ts). \(q\) adalah integer yang mewakili panjang lag yang terbatas.

Pemodelan (Lag=2)

model.dlm <- dlm(x = train$CO,y = train$AQI , q = 2)
summary(model.dlm)

## 
## Call:
## lm(formula = model.formula, data = design)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.3193 -0.7306  0.3403  0.9480  3.5824 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -50.53801   14.21921  -3.554 0.000828 ***
## x.t           0.44406    0.16443   2.701 0.009370 ** 
## x.1          -0.03462    0.26662  -0.130 0.897193    
## x.2          -0.02433    0.16354  -0.149 0.882324    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.634 on 51 degrees of freedom
## Multiple R-squared:  0.4581, Adjusted R-squared:  0.4262 
## F-statistic: 14.37 on 3 and 51 DF,  p-value: 6.574e-07
## 
## AIC and BIC values for the model:
##        AIC      BIC
## 1 215.9118 225.9485

AIC(model.dlm)

## [1] 215.9118

BIC(model.dlm)

## [1] 225.9485

Dari hasil diatas, didapat bahwa \(P-value\) dari intercept (0.000828) dan \(x_{t} (0.009370)<0.05\). Hal ini menunjukkan bahwa intercept dan \(x_{t}\) berpengaruh signifikan terhadap \(y\). Adapun model keseluruhan yang terbentuk adalah sebagai berikut

\[ \hat{Y_t}=-50.53801+0.44406X_t-0.03462X_{t-1}-0.02433X_{t-2} \]

Peramalan dan Akurasi

Berikut merupakan hasil peramalan \(y\) untuk 15 periode kedepan

fore.dlm <- forecast(model = model.dlm, x=test$CO, h=15)
fore.dlm

## $forecasts
##  [1] 25.94516 25.20406 23.77964 23.93581 23.27592 23.33370 23.37430 23.37430
##  [9] 23.37430 23.37430 24.11540 24.79872 25.44144 26.08415 25.98577
## 
## $call
## forecast.dlm(model = model.dlm, x = test$CO, h = 15)
## 
## attr(,"class")
## [1] "forecast.dlm" "dLagM"

mape.dlm <- MAPE(fore.dlm$forecasts, test$AQI)
mape.dlm

## [1] 0.1213772

#akurasi data training
GoF(model.dlm)

##            n      MAE          MPE       MAPE      sMAPE     MASE      MSE
## model.dlm 55 1.158733 -0.003304009 0.04292434 0.04229426 1.862249 2.474342
##                 MRAE    GMRAE
## model.dlm 5069626475 29978.66

Pada perhitungan keakuratan model menggunakan metode Regression with Distributed Lag didapatkan nilai MAPE 4.23%. Nilai akurasi model ini kurang dari 10% sehingga dapat dikategorikan sangat baik.

Lag Optimum

#penentuan lag optimum 
finiteDLMauto(formula = AQI ~ CO,
              data = data.frame(train), q.min = 1, q.max = 6,
              model.type = "dlm", error.type = "AIC", trace = FALSE)

##   q - k    MASE      AIC      BIC    GMRAE    MBRAE R.Adj.Sq    Ljung-Box
## 6     6 1.44155 190.2107 207.5971 18987.97 -0.61847  0.58203 2.058924e-05

Berdasarkan output tersebut, lag optimum didapatkan ketika lag=6. Selanjutnya dilakukan pemodelan untuk lag=6

#model dlm dengan lag optimum
model.dlm2 <- dlm(x = train$CO,y = train$AQI , q = 6)
summary(model.dlm2)

## 
## Call:
## lm(formula = model.formula, data = design)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.7262 -0.4503  0.1114  0.8057  3.5954 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -49.53042   15.68004  -3.159   0.0029 ** 
## x.t           0.75821    0.17159   4.419 6.62e-05 ***
## x.1          -0.17542    0.28067  -0.625   0.5353    
## x.2          -0.30631    0.27110  -1.130   0.2648    
## x.3           0.14117    0.29204   0.483   0.6313    
## x.4           0.22282    0.27520   0.810   0.4226    
## x.5          -0.27209    0.28147  -0.967   0.3391    
## x.6           0.01289    0.16324   0.079   0.9374    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.426 on 43 degrees of freedom
## Multiple R-squared:  0.6405, Adjusted R-squared:  0.582 
## F-statistic: 10.95 on 7 and 43 DF,  p-value: 7.656e-08
## 
## AIC and BIC values for the model:
##        AIC      BIC
## 1 190.2107 207.5971

AIC(model.dlm2)

## [1] 190.2107

BIC(model.dlm2)

## [1] 207.5971

Dari hasil tersebut terdapat peubah yang berpengaruh signifikan terhadap taraf nyata 5% yaitu \(x_t\). Adapun keseluruhan model yang terbentuk adalah

\[ \hat{Y_t}=-49.53042+0.75821X_t-0.17542X_{t-1}-0.30631X_{t-2}+0.14117X_{t-3}+0.22282X_{t-4}-0.27209X_{t-5}+0.01289X_{t-6} \]

Adapun hasil peramalan 15 periode kedepan menggunakan model tersebut adalah sebagai berikut

#peramalan dan akurasi
fore.dlm2 <- forecast(model = model.dlm2, x=test$CO, h=15) #ramal 15 periode kedepan
fore.dlm2

## $forecasts
##  [1] 26.13086 24.49358 22.70963 23.78485 23.30628 22.75598 22.97753 23.62861
##  [9] 23.21371 23.66781 24.91170 25.88434 26.34577 27.04279 26.84629
## 
## $call
## forecast.dlm(model = model.dlm2, x = test$CO, h = 15)
## 
## attr(,"class")
## [1] "forecast.dlm" "dLagM"

mape.dlm2<- MAPE(fore.dlm2$forecasts, test$AQI)
#akurasi data training
GoF(model.dlm2)

##             n       MAE          MPE       MAPE      sMAPE     MASE      MSE
## model.dlm2 51 0.8937606 -0.002271097 0.03307414 0.03264319 1.441549 1.713954
##                  MRAE    GMRAE
## model.dlm2 3249579156 18987.97

Model tersebut merupakan model yang sangat baik dengan nilai MAPE yang kurang dari 10%, yaitu sebesar 3.3%.

Model Autoregressive

Peubah dependen dipengaruhi oleh peubah independen pada waktu sekarang, serta dipengaruhi juga oleh peubah dependen itu sendiri pada satu waktu yang lalu maka model tersebut disebut autoregressive (Gujarati 2004).

Pemodelan

Pemodelan Autoregressive dilakukan menggunakan fungsi dLagM::ardlDlm() . Fungsi tersebut akan menerapkan autoregressive berordo \((p,q)\) dengan satu prediktor. Fungsi umum dari ardlDlm() adalah sebagai berikut.

ardlDlm(formula = NULL , data = NULL , x = NULL , y = NULL , p = 1 , q = 1 , 
         remove = NULL )

Dengan \(p\) adalah integer yang mewakili panjang lag yang terbatas dan \(q\) adalah integer yang merepresentasikan ordo dari proses autoregressive.

model.ardl <- ardlDlm(x = train$CO, y = train$AQI, p = 1 , q = 1)
summary(model.ardl)

## 
## Time series regression with "ts" data:
## Start = 2, End = 57
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.2430 -0.4403  0.1177  0.4929  1.8901 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -5.40310    7.42525  -0.728  0.47008    
## X.t          0.23268    0.06841   3.401  0.00130 ** 
## X.1         -0.18465    0.06709  -2.753  0.00812 ** 
## Y.1          0.83951    0.06939  12.098  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8499 on 52 degrees of freedom
## Multiple R-squared:  0.8511, Adjusted R-squared:  0.8425 
## F-statistic: 99.04 on 3 and 52 DF,  p-value: < 2.2e-16

AIC(model.ardl)

## [1] 146.562

BIC(model.ardl)

## [1] 156.6888

Dari hasil tersebut, didapat bahwa peubah \(x_t (0.00130)\) , \(x_{t-1}(0.00812)\) , dan \(y_{t-1}(< 2e-16)\) memiliki nilai P-Value < 0.05 Hal ini menunjukkan bahwa ketiga peubah tersebut berpengaruh signifikan terhadap \(y_t\) pada taraf nyata 5%. Model keseluruhannya adalah sebagai berikut:

\[ \hat{Y}=-5.40310+0.23268X_t-0.18465X_{t-1}+0.83951Y_{t-1} \]

Peramalan dan Akurasi

fore.ardl <- forecast(model = model.ardl, x=test$CO, h=15)
fore.ardl

## $forecasts
##  [1] 25.12311 24.83814 24.13041 24.15262 23.78294 23.78076 23.77893 23.77739
##  [9] 23.77611 23.77502 24.16245 24.56785 24.98834 25.42150 25.47697
## 
## $call
## forecast.ardlDlm(model = model.ardl, x = test$CO, h = 15)
## 
## attr(,"class")
## [1] "forecast.ardlDlm" "dLagM"

Data di atas merupakan hasil peramalan untuk 15 periode ke depan menggunakan Model Autoregressive dengan \(p=1\) dan \(q=1\).

mape.ardl <- MAPE(fore.ardl$forecasts, test$AQI)
mape.ardl

## [1] 0.1317563

#akurasi data training
GoF(model.ardl)

##             n       MAE           MPE       MAPE      sMAPE      MASE       MSE
## model.ardl 56 0.6377206 -0.0008738902 0.02353982 0.02347038 0.9964384 0.6708097
##                  MRAE    GMRAE
## model.ardl 1418962508 13514.21

Berdasarkan akurasi di atas, terlihat bahwa nilai MAPE keduanya tidak jauh berbeda. Artinya, model regresi dengan distribusi lag ini tidak overfitted atau underfitted

Lag Optimum

#penentuan lag optimum
model.ardl.opt <- ardlBoundOrders(data = data.frame(data), ic = "AIC", 
                                  formula = AQI ~ CO )
min_p=c()
for(i in 1:6){
  min_p[i]=min(model.ardl.opt$Stat.table[[i]])
}
q_opt=which(min_p==min(min_p, na.rm = TRUE))
p_opt=which(model.ardl.opt$Stat.table[[q_opt]] == 
              min(model.ardl.opt$Stat.table[[q_opt]], na.rm = TRUE))
data.frame("q_optimum" = q_opt, "p_optimum" = p_opt, 
           "AIC"=model.ardl.opt$min.Stat)

##   q_optimum p_optimum      AIC
## 1         6        15 121.7055

Dari tabel di atas, dapat terlihat bahwa nilai AIC terendah didapat ketika \(p=15\) dan \(q=6\), yaitu sebesar 121.7055. Artinya, model autoregressive optimum didapat ketika \(p=15\) dan \(q=6\).

Selanjutnya dapat dilakukan pemodelan dengan nilai \(p\) dan \(q\) optimum seperti inisialisasi di langkah sebelumnya.

Pemodelan DLM & ARDL dengan Library `dynlm`

Pemodelan regresi dengan peubah lag tidak hanya dapat dilakukan dengan fungsi pada packages dLagM , tetapi terdapat packages dynlm yang dapat digunakan. Fungsi dynlm secara umum adalah sebagai berikut.

dynlm(formula, data, subset, weights, na.action, method = "qr",
  model = TRUE, x = FALSE, y = FALSE, qr = TRUE, singular.ok = TRUE,
  contrasts = NULL, offset, start = NULL, end = NULL, ...)

Untuk menentukan formula model yang akan digunakan, tersedia fungsi tambahan yang memungkinkan spesifikasi dinamika (melalui d() dan L()) atau pola linier/siklus dengan mudah (melalui trend(), season(), dan harmon()). Semua fungsi formula baru mengharuskan argumennya berupa objek deret waktu (yaitu, "ts" atau "zoo").

#sama dengan model dlm q=1
cons_lm1 <- dynlm(AQI ~ CO+L(CO),data = train.ts)
#sama dengan model ardl p=1 q=0
cons_lm2 <- dynlm(AQI ~ CO+L(AQI),data = train.ts)
#sama dengan ardl p=1 q=1
cons_lm3 <- dynlm(AQI ~ CO+L(CO)+L(AQI),data = train.ts)
#sama dengan dlm p=2
cons_lm4 <- dynlm(AQI ~ CO+L(CO)+L(CO,2),data = train.ts)

Ringkasan Model

summary(cons_lm1)

## 
## Time series regression with "ts" data:
## Start = 2, End = 57
## 
## Call:
## dynlm(formula = AQI ~ CO + L(CO), data = train.ts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.3285 -0.7628  0.3010  0.9793  3.5702 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -48.66712   12.58867  -3.866 0.000304 ***
## CO            0.43680    0.12826   3.406 0.001266 ** 
## L(CO)        -0.06074    0.12826  -0.474 0.637772    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.644 on 53 degrees of freedom
## Multiple R-squared:  0.4319, Adjusted R-squared:  0.4104 
## F-statistic: 20.14 on 2 and 53 DF,  p-value: 3.11e-07

summary(cons_lm2)

## 
## Time series regression with "ts" data:
## Start = 2, End = 57
## 
## Call:
## dynlm(formula = AQI ~ CO + L(AQI), data = train.ts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.4394 -0.6878  0.1422  0.4470  2.1296 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -11.21330    7.54764  -1.486    0.143    
## CO            0.08079    0.04287   1.884    0.065 .  
## L(AQI)        0.81034    0.07271  11.145 1.66e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9011 on 53 degrees of freedom
## Multiple R-squared:  0.8294, Adjusted R-squared:  0.8229 
## F-statistic: 128.8 on 2 and 53 DF,  p-value: < 2.2e-16

summary(cons_lm3)

## 
## Time series regression with "ts" data:
## Start = 2, End = 57
## 
## Call:
## dynlm(formula = AQI ~ CO + L(CO) + L(AQI), data = train.ts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.2430 -0.4403  0.1177  0.4929  1.8901 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -5.40310    7.42525  -0.728  0.47008    
## CO           0.23268    0.06841   3.401  0.00130 ** 
## L(CO)       -0.18465    0.06709  -2.753  0.00812 ** 
## L(AQI)       0.83951    0.06939  12.098  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8499 on 52 degrees of freedom
## Multiple R-squared:  0.8511, Adjusted R-squared:  0.8425 
## F-statistic: 99.04 on 3 and 52 DF,  p-value: < 2.2e-16

summary(cons_lm4)

## 
## Time series regression with "ts" data:
## Start = 3, End = 57
## 
## Call:
## dynlm(formula = AQI ~ CO + L(CO) + L(CO, 2), data = train.ts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.3193 -0.7306  0.3403  0.9480  3.5824 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -50.53801   14.21921  -3.554 0.000828 ***
## CO            0.44406    0.16443   2.701 0.009370 ** 
## L(CO)        -0.03462    0.26662  -0.130 0.897193    
## L(CO, 2)     -0.02433    0.16354  -0.149 0.882324    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.634 on 51 degrees of freedom
## Multiple R-squared:  0.4581, Adjusted R-squared:  0.4262 
## F-statistic: 14.37 on 3 and 51 DF,  p-value: 6.574e-07