1
load("USA_data.RData")
attach(USA_data)
yo<-ts(YO, start= 1960, frequency=1)
io<-ts(IO, start= 1960, frequency=1)
lnyo<-log(yo)
lnio<-log(io)
llnio <- embed(lnio, dimension = 2)
head(llnio)
## [,1] [,2]
## [1,] 13.34799 13.31652
## [2,] 13.44161 13.34799
## [3,] 13.49743 13.44161
## [4,] 13.56089 13.49743
## [5,] 13.65956 13.56089
## [6,] 13.74639 13.65956
llnio0 <- llnio[, 1]
llnio1 <- llnio[, 2]
llnyo<-embed(lnyo, dimension =2)
llnyo0<-llnyo[, 1]
llnyo1<-llnyo[, 2]
2
If two series same stochastic trend and move together in the long
run, they are cointegrated.
3
#a
library(bootUR)
LNAPIO <- lnio - lnyo
# 1, just sending lnyo to other side.
#c
ts.plot(LNAPIO, ylab = "LNAPIO")
#Not stationary
#d
adfLNAPIO <- adf(LNAPIO, deterministics = "trend")
adfLNAPIO
##
## Two-step ADF test (with trend) on a single time series
##
## data: LNAPIO
## null hypothesis: Series has a unit root
## alternative hypothesis: Series is stationary
##
## estimate largest root statistic p-value
## LNAPIO 0.7198 -3.198 0.09617
#I included trend because the series seemed to trend upwards.
4
#a
#If it is spurious or not
#b
r19 <- lm(lnio ~ lnyo)
summary(r19)
##
## Call:
## lm(formula = lnio ~ lnyo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.17891 -0.03535 0.01693 0.04482 0.09827
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.31213 0.23770 -18.14 <2e-16 ***
## lnyo 1.16058 0.01465 79.21 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06309 on 62 degrees of freedom
## Multiple R-squared: 0.9902, Adjusted R-squared: 0.9901
## F-statistic: 6274 on 1 and 62 DF, p-value: < 2.2e-16
#c
res19 <- residuals(r19)
ts.plot(res19, ylab = "Residuals lnIO ~ lnYO")
# Maybe a weak statinonary, but not really
#d
adfresid <- adf(res19, deterministics = "none")
adfresid
##
## Two-step ADF test (with none) on a single time series
##
## data: res19
## null hypothesis: Series has a unit root
## alternative hypothesis: Series is stationary
##
## estimate largest root statistic p-value
## res19 0.7427 -3.041 0.003027
#e
#resiudals are estimated and the nature of MLR's force it to demonstrate stationarity properties like zeroo conditional mean.
#f
#less negative cannot reject null
#g residuals are estimated, thus of course extra noise is added. Q3 test is more powrful.
5
r20 <- lm(llnio0~llnyo0+llnyo1+llnio1)
summary(r20)
##
## Call:
## lm(formula = llnio0 ~ llnyo0 + llnyo1 + llnio1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.071637 -0.012508 0.003956 0.012111 0.080946
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.5464 0.2457 -6.294 4.18e-08 ***
## llnyo0 2.9219 0.1681 17.382 < 2e-16 ***
## llnyo1 -2.5986 0.1764 -14.731 < 2e-16 ***
## llnio1 0.7416 0.0517 14.342 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02561 on 59 degrees of freedom
## Multiple R-squared: 0.9984, Adjusted R-squared: 0.9983
## F-statistic: 1.207e+04 on 3 and 59 DF, p-value: < 2.2e-16
# impact hat = 2.9219, impact = 2.9219, they are same no need to plug back, sd= 0.1681
B <- coef(r20)
mr20<-(B["llnyo0"]+ B["llnyo1"])+(B["llnio1"]*B["llnyo0"])
mr20
## llnyo0
## 2.490196
nls_theta2 = nls(llnio0 ~ beta0 + ((theta2 -beta2)/(1+ beta3))*llnyo0 + beta2*llnyo1 +
beta3*llnio1, start = list(beta0 = 1, theta2 = 1, beta2 = 1, beta3 = 1))
summary(nls_theta2)
##
## Formula: llnio0 ~ beta0 + ((theta2 - beta2)/(1 + beta3)) * llnyo0 + beta2 *
## llnyo1 + beta3 * llnio1
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## beta0 -1.5464 0.2457 -6.294 4.18e-08 ***
## theta2 2.4902 0.1562 15.948 < 2e-16 ***
## beta2 -2.5986 0.1764 -14.731 < 2e-16 ***
## beta3 0.7416 0.0517 14.342 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02561 on 59 degrees of freedom
##
## Number of iterations to convergence: 2
## Achieved convergence tolerance: 4.361e-06
# medium= 2.4902, sd = 0.1562
lr20<- (B["llnyo0"]+ B["llnyo1"])/(1-B["llnio1"])
lr20 #long = 1.251296
## llnyo0
## 1.251296
nls_theta3 = nls(llnio0 ~ beta0 + (theta3-theta3*beta3-beta2)*llnyo0 + beta2*llnyo1 +
beta3*llnio1, start = list(beta0 = -1.5464, theta3 = 1.251296, beta2 = -2.5986, beta3 = 0.7416 ))
summary(nls_theta3)
##
## Formula: llnio0 ~ beta0 + (theta3 - theta3 * beta3 - beta2) * llnyo0 +
## beta2 * llnyo1 + beta3 * llnio1
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## beta0 -1.54642 0.24568 -6.294 4.18e-08 ***
## theta3 1.25130 0.03102 40.342 < 2e-16 ***
## beta2 -2.59857 0.17640 -14.731 < 2e-16 ***
## beta3 0.74157 0.05170 14.342 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02561 on 59 degrees of freedom
##
## Number of iterations to convergence: 1
## Achieved convergence tolerance: 1.412e-06
#long= 1.25130, sd = 0.03102
6
r21 <- lm(llnio0~llnyo0+llnyo1)
summary(r21)
##
## Call:
## lm(formula = llnio0 ~ llnyo0 + llnyo1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.11652 -0.03579 0.01002 0.04068 0.11918
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.7245 0.2229 -21.199 < 2e-16 ***
## llnyo0 2.9272 0.3531 8.290 1.57e-11 ***
## llnyo1 -1.7444 0.3488 -5.002 5.25e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05378 on 60 degrees of freedom
## Multiple R-squared: 0.9927, Adjusted R-squared: 0.9925
## F-statistic: 4081 on 2 and 60 DF, p-value: < 2.2e-16
# impact hat = 2.9219, impact = 2.9219, they are same no need to plug back, sd= 0.1681
C <- coef(r21)
mr21<-(C["llnyo0"]+ C["llnyo1"])
mr21# = 1.182772
## llnyo0
## 1.182772
nls_theta2 = nls(llnio0 ~ beta0 + (theta2 - beta2)*llnyo0 + beta2*llnyo1,
start = list(beta0 = 1, theta2 = 1, beta2 = 1))
summary(nls_theta2)
##
## Formula: llnio0 ~ beta0 + (theta2 - beta2) * llnyo0 + beta2 * llnyo1
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## beta0 -4.72451 0.22287 -21.199 < 2e-16 ***
## theta2 1.18277 0.01352 87.491 < 2e-16 ***
## beta2 -1.74443 0.34876 -5.002 5.25e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05378 on 60 degrees of freedom
##
## Number of iterations to convergence: 1
## Achieved convergence tolerance: 1.421e-06
# medium= 1.18277, sd = 0.01352
lr21<- (C["llnyo0"]+ C["llnyo1"])
lr21
## llnyo0
## 1.182772
#long = 1.251296
nls_theta3 = nls(llnio0 ~ beta0 + (theta3-beta2)*llnyo0 + beta2*llnyo1
, start = list(beta0 = 1, theta3 =1, beta2 = 1))
summary(nls_theta3)
##
## Formula: llnio0 ~ beta0 + (theta3 - beta2) * llnyo0 + beta2 * llnyo1
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## beta0 -4.72451 0.22287 -21.199 < 2e-16 ***
## theta3 1.18277 0.01352 87.491 < 2e-16 ***
## beta2 -1.74443 0.34876 -5.002 5.25e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05378 on 60 degrees of freedom
##
## Number of iterations to convergence: 1
## Achieved convergence tolerance: 1.421e-06
#long= 1.18277, sd = 0.01352
7
#a
# We can repeadealty substitute y lag which has coefficent smaller than one, thus this lead to a dynamic model.
r22 <- lm(llnio0~llnyo0+llnio1)
summary(r22)
##
## Call:
## lm(formula = llnio0 ~ llnyo0 + llnio1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.194138 -0.018083 0.008899 0.037521 0.112631
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.1418 0.5198 -4.121 0.000118 ***
## llnyo0 0.5944 0.1231 4.829 9.85e-06 ***
## llnio1 0.4844 0.1044 4.641 1.93e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05492 on 60 degrees of freedom
## Multiple R-squared: 0.9924, Adjusted R-squared: 0.9921
## F-statistic: 3913 on 2 and 60 DF, p-value: < 2.2e-16
# impact hat = 2.9219, impact = 2.9219, they are same no need to plug back, sd= 0.1681
D <- coef(r22)
mr22<-D["llnyo0"]+(D["llnyo0"]* D["llnio1"])
mr22#=0.8822963
## llnyo0
## 0.8822963
nls_theta2 = nls(llnio0 ~ beta0 + (theta2/(1+beta3))*llnyo0 + beta3*llnio1,
start = list(beta0 = 1, theta2 = 1, beta3 = 1))
summary(nls_theta2)
##
## Formula: llnio0 ~ beta0 + (theta2/(1 + beta3)) * llnyo0 + beta3 * llnio1
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## beta0 -2.1418 0.5198 -4.121 0.000118 ***
## theta2 0.8823 0.1212 7.279 8.35e-10 ***
## beta3 0.4844 0.1044 4.641 1.93e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05492 on 60 degrees of freedom
##
## Number of iterations to convergence: 2
## Achieved convergence tolerance: 4.75e-07
# medium= 0.8823, sd = 0.1212
D <- coef(r22)
lr22<- D["llnyo0"]/(1-D["llnio1"])
lr22 #long = 1.152825
## llnyo0
## 1.152825
nls_theta3 = nls(llnio0 ~ beta0 + (theta3-theta3*beta3)*llnyo0 + beta3*llnio1,
start = list(beta0 = -1.5464, theta3 = 1.251296, beta3 = 0.7416 ))
summary(nls_theta3)
##
## Formula: llnio0 ~ beta0 + (theta3 - theta3 * beta3) * llnyo0 + beta3 *
## llnio1
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## beta0 -2.14183 0.51976 -4.121 0.000118 ***
## theta3 1.15282 0.02583 44.639 < 2e-16 ***
## beta3 0.48442 0.10439 4.641 1.93e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05492 on 60 degrees of freedom
##
## Number of iterations to convergence: 7
## Achieved convergence tolerance: 1.22e-06
#long= 1.15282, sd = 0.02583
8
#a
#If they are cointegrated it is stationary.
r23r <- lm(lnio~lnyo)
res<- r23r$residuals
res <- res[1:(length(res) - 1)]
dlnio <- diff(lnio)
dlnyo <- diff(lnyo)
length(dlnio)
## [1] 63
length(res)
## [1] 63
r23 <- lm(dlnio~ dlnyo+res)
summary(r23)
##
## Call:
## lm(formula = dlnio ~ dlnyo + res)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.087454 -0.018205 0.001911 0.017386 0.087842
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.046115 0.006251 -7.378 5.67e-10 ***
## dlnyo 2.715524 0.173684 15.635 < 2e-16 ***
## res -0.260980 0.056739 -4.600 2.24e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0281 on 60 degrees of freedom
## Multiple R-squared: 0.8152, Adjusted R-squared: 0.809
## F-statistic: 132.3 on 2 and 60 DF, p-value: < 2.2e-16
9
ecm1 <- lnio[1:(length(lnio)-1)] - lnyo[1:(length(lnyo)-1)]
r24 <- lm(dlnio ~ dlnyo + ecm1)
summary(r24)
##
## Call:
## lm(formula = dlnio ~ dlnyo + ecm1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.081970 -0.020743 0.001277 0.020289 0.075018
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.054401 0.068343 -0.796 0.429
## dlnyo 2.702614 0.209762 12.884 <2e-16 ***
## ecm1 -0.005158 0.040694 -0.127 0.900
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03268 on 60 degrees of freedom
## Multiple R-squared: 0.75, Adjusted R-squared: 0.7417
## F-statistic: 90.02 on 2 and 60 DF, p-value: < 2.2e-16
alpha0 <- coef(r24)["(Intercept)"]
alpha1 <- coef(r24)["dlnyo"]
delta <- coef(r24)["ecm1"]
alpha0
## (Intercept)
## -0.05440105
alpha1
## dlnyo
## 2.702614
delta
## ecm1
## -0.005158375
10
r25 <- lm(dlnio ~ dlnyo)
summary(r25)
##
## Call:
## lm(formula = dlnio ~ dlnyo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.081214 -0.021075 0.001275 0.019652 0.074515
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.045787 0.007209 -6.351 2.98e-08 ***
## dlnyo 2.709796 0.200330 13.527 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03241 on 61 degrees of freedom
## Multiple R-squared: 0.75, Adjusted R-squared: 0.7459
## F-statistic: 183 on 1 and 61 DF, p-value: < 2.2e-16
impact25 <- coef(r25)[2]
medium25 <- impact25
11
impact26 <- coef(r19)[2]
medium26 <- impact26
long26 <- impact26
12
APIO <- IO / YO
APIOs <- exp(-alpha0 / delta)
meanAPIO <- mean(APIO)
rangeAPIO <- range(APIO)
APIOs
## (Intercept)
## 2.629423e-05
meanAPIO
## [1] 0.1821689
rangeAPIO
## [1] 0.1431088 0.2180714