Tải dữ liệu
library(tsDyn)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
data(lynx)
head(lynx)
## Time Series:
## Start = 1821
## End = 1826
## Frequency = 1
## [1] 269 321 585 871 1475 2821
Test of linearity against ( TÌm số ngưỡng)
kiemtra1 <- setarTest(lynx, m=1,thDelay=0, test=c("1vs"))
summary(kiemtra1)
## Test of linearity against setar(2) and setar(3)
##
## Test Pval
## 1vs2 20.05733 0
## 1vs3 42.81466 0
##
## Critical values:
## 0.9 0.95 0.975 0.99
## 1vs2 11.53407 12.08120 12.35476 12.51890
## 1vs3 17.03573 17.67661 17.99706 18.18933
##
## SSR of original series:
## SSR
## AR 137160013
## SETAR(2) 116484232
## SETAR(3) 99471266
##
## Threshold of original series:
## th1 th2
## SETAR(2) 3465 NA
## SETAR(3) 1676 3465
##
## Number of bootstrap replications: 10
kiemtra2 <- setarTest(lynx, m=1,thDelay=0, test=c("2vs3"))
## Warning: Possible unit root in the regime. Roots are: 0.9425
## Warning: Possible unit root in the regime. Roots are: 0.9443
## Warning: Possible unit root in the regime. Roots are: 0.9425
summary(kiemtra2)
## Test of setar(2) against setar(3)
##
## Test Pval
## 2vs3 19.32684 0
##
## Critical values:
## 0.9 0.95 0.975 0.99
## 2vs3 9.992291 10.79291 11.19322 11.43341
##
## SSR of original series:
## SSR
## AR 137160013
## SETAR(2) 116484232
## SETAR(3) 99471266
##
## Threshold of original series:
## th1 th2
## SETAR(2) 3465 NA
## SETAR(3) 1676 3465
##
## Number of bootstrap replications: 10
Tìm giá trị ngưỡng
grid <- selectSETAR(lynx, m=1, thDelay=0,trim= 0.15, criterion="SSR", nthresh=2)
## Using maximum autoregressive order for low regime: mL = 1
## Using maximum autoregressive order for high regime: mH = 1
## Using maximum autoregressive order for middle regime: mM = 1
## Searching on 75 possible threshold values within regimes with sufficient ( 15% ) number of observations
## Searching on 75 combinations of thresholds (75) and thDelay (1)
## Result of the one threshold search:
## -Thresh: 1388 -Delay: 0 - SSR 123102676
## Second best: 2577 (conditionnal on th= 1388 and Delay= 0 ) SSR/AIC: 114452658
## Second best: 1000 (conditionnal on th= 2577 and Delay= 0 ) SSR/AIC: 113310032
summary(grid)
## Length Class Mode
## res 3 data.frame list
## res2 3 data.frame list
## bests 4 -none- numeric
## th 2 -none- numeric
## nthresh 1 -none- numeric
## allTh 225 -none- numeric
## criterion 1 -none- character
## firstBests 3 -none- numeric
## bests2th 3 -none- numeric
## ML 1 -none- numeric
## MM 1 -none- numeric
## MH 1 -none- numeric
## same.lags 1 -none- logical
grid
## Results of the grid search for 1 threshold
## Conditional on m= 1
## thDelay th SSR
## 1 0 1388 123102676
## 2 0 1307 123951941
## 3 0 1475 124388924
## 4 0 1676 124516444
## 5 0 1638 124557228
##
## Results of the grid search for 2 thresholds
## Conditional on thDelay = 0 and m = 1
## th1 th2 SSR
## 1 1000 2577 113310032
##
## Overall best results:
## thDelay th1 th2 SSR
## 0 1000 2577 113310032
## With lags:
## -ML: 1
## -MM: 1
## -MH: 1
Ước lượng ngưỡng
set <-setar(lynx, m=1,thDelay=0, th=grid$th)
## Warning: Possible unit root in the low regime. Roots are: 0.7924
summary(set)
##
## Non linear autoregressive model
##
## SETAR model ( 3 regimes)
## Coefficients:
## Low regime:
## const.L phiL.1
## 84.715767 1.261983
##
## Mid regime:
## const.M phiM.1
## 3612.9339667 -0.8987912
##
## High regime:
## const.H phiH.1
## 2098.6793167 0.3328867
##
## Threshold:
## -Variable: Z(t) = + (1) X(t)
## -Value: 1000 2577 (fixed)
## Proportion of points in low regime: 54.87% Middle regime: 21.24% High regime: 23.89%
##
## Residuals:
## Min 1Q Median 3Q Max
## -2827.779 -376.600 -27.809 347.588 3738.868
##
## Fit:
## residuals variance = 993948, AIC = 1586, MAPE = 105.3%
##
## Coefficient(s):
##
## Estimate Std. Error t value Pr(>|t|)
## const.L 84.71577 241.77064 0.3504 0.726723
## phiL.1 1.26198 0.52230 2.4162 0.017361 *
## const.M 3612.93397 856.55362 4.2180 5.144e-05 ***
## phiM.1 -0.89879 0.46090 -1.9501 0.053760 .
## const.H 2098.67932 644.49843 3.2563 0.001508 **
## phiH.1 0.33289 0.15855 2.0996 0.038097 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold
## Variable: Z(t) = + (1) X(t)
##
## Value: 1000 2577 (fixed)
Chạy tìm số ngưỡng và giá trị TỰ ĐỘNG
mod.star <- star(log10(lynx), mTh=c(0,1), control=list(maxit=3000))
## Testing linearity... p-Value = 0.0001831653
## The series is nonlinear. Incremental building procedure:
## Building a 2 regime STAR.
## Performing grid search for starting values...
## Starting values fixed: gamma = 11.15385 , th = 3.337486 ; SSE = 4.337664
## Optimization algorithm converged
## Optimized values fixed for regime 2 : gamma = 11.15383 , th = 3.339199 ; SSE = 4.337643
##
## Testing for addition of regime 3.
## Estimating gradient matrix...
## Done. Computing the test statistic...
## Done. Regime 3 is NOT accepted (p-Value = 0.9843362 ).
##
## Finished building a MRSTAR with 2 regimes
mod.star
##
## Non linear autoregressive model
##
## Multiple regime STAR model
##
## Regime 1 :
## Linear parameters: 0.4891014, 1.2465399, -0.3664328
##
## Regime 2 :
## Linear parameters: -1.0240758, 0.4232669, -0.2546088
## Non-linear parameters:
## 11.1538344, 3.3391985
summary(mod.star)
##
## Non linear autoregressive model
##
## Multiple regime STAR model
##
## Regime 1 :
## Linear parameters: 0.4891014, 1.2465399, -0.3664328
##
## Regime 2 :
## Linear parameters: -1.0240758, 0.4232669, -0.2546088
## Non-linear parameters:
## 11.1538344, 3.3391985
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.59482018 -0.10735988 0.01430922 0.11109839 0.51034240
##
## Fit:
## residuals variance = 0.03805, AIC = -357, MAPE = 5.58%
addRegime(mod.star)
##
## Non linear autoregressive model
##
## Multiple regime STAR model
##
## Regime 1 :
## Linear parameters: 0.4891014, 1.2465399, -0.3664328
##
## Regime 2 :
## Linear parameters: -1.0240758, 0.4232669, -0.2546088
## Non-linear parameters:
## 11.1538344, 3.3391985
## Regime 3 :
## Linear parameters: -0.9388364, 1.3367467, -0.0239093
## Non-linear parameters:
## 0, 0
Multivariate Threshold Vector Autoregressive model
data(zeroyld)
head(zeroyld)
## short.run long.run
## 1 2.183 1.575
## 2 2.246 1.545
## 3 2.308 1.762
## 4 2.401 1.773
## 5 2.558 1.808
## 6 2.507 1.810
tv <- TVAR(zeroyld, lag=2, nthresh=2, thDelay=1, trim=0.1, mTh=1, plot=FALSE)
## Best unique threshold 10.653
## Second best: 8.124 (conditionnal on th= 10.653 and Delay= 1 ) SSR/AIC: 164.1456
## Second best: 10.653 (conditionnal on th= 8.124 and Delay= 1 ) SSR/AIC: 164.1456
##
## Second step best thresholds 8.124 10.653 SSR: 164.1456
print(tv)
## Model TVAR with 2 thresholds
##
## $Bdown
## Intercept short.run -1 long.run -1 short.run -2
## Equation short.run 0.01496509 1.0061051 0.03855754 -0.007405335
## Equation long.run -0.03467035 0.3622364 1.02090855 -0.287156749
## long.run -2
## Equation short.run -0.03455587
## Equation long.run -0.09201102
##
## $Bmiddle
## Intercept short.run -1 long.run -1 short.run -2 long.run -2
## Equation short.run -0.2169119 0.83059144 0.0127026 0.13505057 0.04939883
## Equation long.run -0.1386902 -0.02142815 0.8490303 0.02085644 0.16932567
##
## $Bup
## Intercept short.run -1 long.run -1 short.run -2 long.run -2
## Equation short.run 1.561918 1.1780904 -0.03076567 -0.3468208 0.06986963
## Equation long.run 2.037192 0.7989221 0.81823826 -0.8562527 0.05771228
##
##
## Threshold value[1] "8.124 10.653"
summary(tv)
## Model TVAR with 2 thresholds
##
## Full sample size: 482 End sample size: 480
## Number of variables: 2 Number of estimated parameters: 30 + 2
## AIC -2155.439 BIC -2021.878 SSR 164.1456
##
## [[1]]
## Intercept short.run -1 long.run -1
## Equation short.run 0.0150(0.0491) 1.0061(0.0990)*** 0.0386(0.0595)
## Equation long.run -0.0347(0.0836) 0.3622(0.1685)* 1.0209(0.1014)***
## short.run -2 long.run -2
## Equation short.run -0.0074(0.0969) -0.0346(0.0589)
## Equation long.run -0.2872(0.1650). -0.0920(0.1003)
##
## [[2]]
## Intercept short.run -1 long.run -1
## Equation short.run -0.2169(0.4479) 0.8306(0.1566)*** 0.0127(0.1005)
## Equation long.run -0.1387(0.7625) -0.0214(0.2666) 0.8490(0.1711)***
## short.run -2 long.run -2
## Equation short.run 0.1351(0.1459) 0.0494(0.1013)
## Equation long.run 0.0209(0.2483) 0.1693(0.1724)
##
## [[3]]
## Intercept short.run -1 long.run -1
## Equation short.run 1.5619(0.5115)** 1.1781(0.1349)*** -0.0308(0.0657)
## Equation long.run 2.0372(0.8707)* 0.7989(0.2297)*** 0.8182(0.1118)***
## short.run -2 long.run -2
## Equation short.run -0.3468(0.1249)** 0.0699(0.0640)
## Equation long.run -0.8563(0.2127)*** 0.0577(0.1089)
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold value: 8.124 10.653
## Percentage of Observations in each regime: 71.5% 16.5% 12.1%
# Test
data<-zeroyld
TVAR.LRtest(data, lag=2, mTh=1,thDelay=1:2, nboot=3, plot=FALSE, trim=0.1, test="1vs")
## Warning: the thDelay values do not correspond to the univariate implementation in tsdyn
## $bestDelay
## Var1
## 1
##
## $LRtest.val
## 1vs2 1vs3
## 30.27935 42.97058
##
## $Pvalueboot
## 1vs2 1vs3
## 0.0000000 0.3333333
##
## $CriticalValBoot
## 90% 95% 97.5% 99%
## 1vs2 22.44646 23.07519 23.38956 23.57818
## 1vs3 45.30365 46.27353 46.75847 47.04944
##
## $type
## [1] "1vs"
##
## attr(,"class")
## [1] "TVARtest"
Threshold Vector Error Correction model
tvec <- TVECM(zeroyld, nthresh=2,lag=1, ngridBeta=20, ngridTh=30, plot=TRUE,trim=0.05, common="All")
## 29 (4.8%) points of the grid lead to regimes with percentage of observations < trim and were not computed
## Best threshold from first search -1.312
## Best cointegrating value 1.045747
## Second best (conditionnal on the first one) -1.312 0.774 SSR 155.4941
## Second step best thresholds -1.312 0.774 SSR 155.4941
print(tvec)
## Model TVECM with 2 thresholds
##
## $Bdown
## ECT Const short.run t -1 long.run t -1
## Equation short.run 0.2343191 0.5627267 0.2133872 -0.13045497
## Equation long.run 1.1492230 1.9445867 0.7652367 0.02490849
##
## $Bmiddle
## ECT Const short.run t -1 long.run t -1
## Equation short.run -0.01031875 0.007322243 0.09944158 0.06907443
## Equation long.run 0.01991474 -0.005735839 0.32699066 0.23687489
##
## $Bup
## ECT Const short.run t -1 long.run t -1
## Equation short.run -0.04920032 0.06691590 -0.11216037 0.04637838
## Equation long.run 0.10732160 -0.09079342 0.08524092 -0.02873180
##
##
## Threshold value[1] -1.312 0.774
summary(tvec)
## #############
## ###Model TVECM
## #############
## Full sample size: 482 End sample size: 480
## Number of variables: 2 Number of estimated parameters 24
## AIC -2223.223 BIC -2114.704 SSR 155.4941
##
##
## Cointegrating vector: (1, - 1.045747 )
## $Bdown
## ECT Const short.run t -1
## Equation short.run 0.2343(0.0207)* 0.5627(0.0041)** 0.2134(0.2238)
## Equation long.run 1.1492(1.0e-11)*** 1.9446(2.1e-09)*** 0.7652(0.0077)**
## long.run t -1
## Equation short.run -0.1305(0.1683)
## Equation long.run 0.0249(0.8718)
##
## $Bmiddle
## ECT Const short.run t -1
## Equation short.run -0.0103(0.7758) 0.0073(0.6884) 0.0994(0.3462)
## Equation long.run 0.0199(0.7362) -0.0057(0.8473) 0.3270(0.0580).
## long.run t -1
## Equation short.run 0.0691(0.2101)
## Equation long.run 0.2369(0.0086)**
##
## $Bup
## ECT Const short.run t -1
## Equation short.run -0.0492(0.4402) 0.0669(0.4232) -0.1122(0.2911)
## Equation long.run 0.1073(0.3022) -0.0908(0.5053) 0.0852(0.6227)
## long.run t -1
## Equation short.run 0.0464(0.5256)
## Equation long.run -0.0287(0.8095)
##
##
## Threshold
## Values: -1.312 0.774
## Percentage of Observations in each regime 6.7% 56.9% 36.5%
#test1<-TVECM.HStest(zeroyld, lag=1, intercept=TRUE, nboot=10)
#test1
#summary(test1)
#plot(test1)
#TVECM.SeoTest(zeroyld,lag=2, beta=1, trim=0.1,nboot=10, plot=FALSE,check=FALSE)
Chạy phản ứng xung
#data(barry)
#dulieu <-barry
## For VAR
#mod_var <- lineVar(dulieu, lag = 2)
#irf(mod_var, impulse = "dolcan", response = c("dolcan", "cpiUSA", "cpiCAN"), boot = FALSE)
## For VECM
#mod_VECM <- VECM(dulieu, lag = 2, estim="ML", r=2)
#irf(mod_VECM, impulse = "dolcan", response = c("dolcan", "cpiUSA", "cpiCAN"), boot =FALSE)