This function simulates multivariate time series generated by a user-specified Threshold Autoregressive (TAR) model.
Usage
simtar(
n,
k = 2,
ars = ars(),
Intercept = TRUE,
trend = c("none", "linear", "quadratic"),
nseason = NULL,
parms,
delay = 0,
thresholds = NULL,
t.series = NULL,
ex.series = NULL,
dist = c("Gaussian", "Student-t", "Hyperbolic", "Laplace", "Slash",
"Contaminated normal", "Skew-Student-t", "Skew-normal"),
skewness = NULL,
extra = NULL,
setar = NULL,
Verbose = TRUE
)Arguments
- n
A positive integer specifying the length of the simulated output series.
- k
A positive integer specifying the dimension of the multivariate output series.
- ars
A list defining the TAR model structure, composed of four elements: the number of regimes (
nregim), the autoregressive order (p), and the maximum lags of the exogenous (q) and threshold (d) series within each regime. This object can be validated using thears()function.- Intercept
An optional logical indicating whether an intercept term is included in each regime.
- trend
An optional character string specifying the degree of deterministic time trend to be included in each regime. Available options are a linear trend (
"linear"), a quadratic trend ("quadratic"), or no trend ("none"). By default,trendis set to"none".- nseason
An optional integer, greater than or equal to 2, specifying the number of seasonal periods. When provided,
nseason - 1seasonal dummy variables are added to the regressors within each regime. By default,nseasonis set toNULL, thereby indicating that the TAR model has no seasonal effects.- parms
A list with one sublist per regime. Each sublist contains two matrices: the first matrix corresponds to the location parameters, and the second matrix corresponds to the scale parameters of the model.
- delay
An optional non-negative integer specifying the delay parameter of the threshold series. By default,
delayis set to0.- thresholds
A numeric vector of length
nregim-1containing the threshold values, sorted in ascending order.- t.series
A matrix containing the values of the threshold series.
- ex.series
A matrix containing the values of the multivariate exogenous series.
- dist
An optional character string specifying the multivariate distribution chosen to model the noise process. Supported options include Gaussian (
"Gaussian"), Student-\(t\) ("Student-t"), Slash ("Slash"), Symmetric Hyperbolic ("Hyperbolic"), Laplace ("Laplace"), and Contaminated Normal ("Contaminated normal"). By default,distis set to"Gaussian".- skewness
An optional numeric vector specifying the skewness parameters of the noise distribution, when applicable.
- extra
An optional value specifying the extra parameter of the noise distribution, when required.
- setar
An optional positive integer indicating which component of the output series is used as the threshold variable. By default,
setaris set toNULL, indicating that the model is not of SETAR type.- Verbose
An optional logical indicating whether a description of the simulated TAR model should be printed. By default,
Verboseis set toTRUE.
Value
A data.frame containing the simulated multivariate output series, and, if
specified, the threshold series and multivariate exogenous series.
References
Vanegas, L.H. and Calderón, S.A. and Rondón, L.M. (2025) Bayesian estimation of a multivariate tar model when the noise process distribution belongs to the class of gaussian variance mixtures. International Journal of Forecasting.
Examples
# \donttest{
###### Simulation of a trivariate TAR model with two regimes
n <- 2000
k <- 3
myars <- ars(nregim=2,p=c(1,2))
Z <- as.matrix(arima.sim(n=n+max(myars$p),list(ar=c(0.5))))
probs <- sort((0.6 + runif(myars$nregim-1)*0.8)*c(1:(myars$nregim-1))/myars$nregim)
dist <- "Student-t"; extra <- 6
parms <- list()
for(j in 1:myars$nregim){
np <- 1 + myars$p[j]*k
parms[[j]] <- list()
parms[[j]]$location <- c(ifelse(runif(np*k)<=0.5,1,-1)*rbeta(np*k,shape1=4,shape2=16))
parms[[j]]$location <- matrix(parms[[j]]$location,np,k)
parms[[j]]$scale <- rgamma(k,shape=1,scale=1)*diag(k)
}
thresholds <- quantile(Z,probs=probs)
out1 <- simtar(n=n, k=k, ars=myars, parms=parms, thresholds=thresholds,
t.series=Z, dist=dist, extra=extra)
#>
#>
#> Sample size :2000 time points
#>
#> Output Series :Y1 | Y2 | Y3
#>
#> Threshold Series :Z with a delay equal to 0
#>
#> Error Distribution :Student-t(6)
#>
#> Number of regimes :2
#>
#> Deterministics :Intercept
#>
#> Autoregressive orders:1, 2
#>
#>
#> Thresholds
#>
#> Regime 1 (-Inf,-0.528421387520755]
#> Regime 2 (-0.528421387520755,Inf)
#>
#>
#> Regime 1:
#>
#> Autoregressive coefficients
#> Intercept phi_1
#> Y1 0.1390502 | 0.23450807 0.2395519 -0.3314744 |
#> Y2 0.1614532 | -0.11454529 -0.3416702 0.2721627 |
#> Y3 -0.2362589 | -0.07889005 0.1288829 0.4240007 |
#>
#> Scale parameter
#> Y1 Y2 Y3
#> Y1 0.04501187 0.000000 0.000000
#> Y2 0.00000000 1.285029 0.000000
#> Y3 0.00000000 0.000000 2.142067
#>
#> Regime 2:
#>
#> Autoregressive coefficients
#> Intercept phi_1 phi_2
#> Y1 0.1888158 | 0.2465074 -0.1258429 -0.08526607 | -0.1551975 -0.3517518
#> Y2 0.1929726 | -0.2814762 0.1934772 0.31548840 | -0.1821664 0.1951535
#> Y3 -0.2147909 | 0.1963584 0.1076944 -0.23977121 | -0.2678406 -0.1112578
#>
#> Y1 0.2233804 |
#> Y2 0.1784293 |
#> Y3 0.1989515 |
#>
#> Scale parameter
#> Y1 Y2 Y3
#> Y1 0.1434533 0.0000000 0.000000
#> Y2 0.0000000 0.2859605 0.000000
#> Y3 0.0000000 0.0000000 2.419312
str(out1)
#> 'data.frame': 2002 obs. of 5 variables:
#> $ Y1 : num 1.0437 0.7045 0.7587 -0.0456 -0.2312 ...
#> $ Y2 : num -0.3704 1.4453 0.0124 -1.3255 0.094 ...
#> $ Y3 : num -0.6342 0.0884 0.6423 2.1867 0.9079 ...
#> $ Regime: num NA NA 1 1 2 2 1 1 1 1 ...
#> $ Z : num -2.373 -2.558 -0.927 -1.25 -0.234 ...
fit1 <- mtar(~ Y1 + Y2 + Y3 | Z, data=out1, ars=myars, dist=dist,
n.burn=2000, n.sim=3000, n.thin=2)
summary(fit1)
#>
#>
#> Sample size : 1999 time points
#>
#> Output Series : Y1 | Y2 | Y3
#>
#> Threshold Series : Z with a estimated delay equal to 0
#>
#> Error Distribution : Student-t
#>
#> Number of regimes : 2
#>
#> Deterministics : Intercept
#>
#> Autoregressive orders: 1, 2
#>
#>
#>
#> Thresholds (Mean, HDI.Lower, HDI.Upper)
#>
#> Regime 1 (-Inf,-0.52671] (-Inf,-0.52984] (-Inf,-0.52221]
#> Regime 2 (-0.52671,Inf) (-0.52984,Inf) (-0.52221,Inf)
#>
#>
#> Regime1:
#>
#> Autoregressive coefficients
#> Mean 2(1-PD) HDI.Lower HDI.Upper Mean 2(1-PD)
#> (Intercept) 0.13798 1e-05 0.12062 0.15573 | 0.21445 0.00001
#> Y1.lag(1) 0.21257 1e-05 0.18995 0.23709 | -0.14989 0.03067
#> Y2.lag(1) 0.23702 1e-05 0.22433 0.25071 | -0.32540 0.00001
#> Y3.lag(1) -0.32878 1e-05 -0.33841 -0.31967 | 0.27427 0.00001
#> HDI.Lower HDI.Upper Mean 2(1-PD) HDI.Lower HDI.Upper
#> (Intercept) 0.12098 0.31459 | -0.27844 0.00001 -0.38780 -0.16161
#> Y1.lag(1) -0.28375 -0.02143 | -0.14516 0.06133 -0.30251 0.00289
#> Y2.lag(1) -0.39751 -0.25560 | 0.23447 0.00001 0.15264 0.32048
#> Y3.lag(1) 0.21938 0.32450 | 0.40656 0.00001 0.34342 0.47198
#>
#> Scale parameter (Mean, HDI.Lower, HDI.Upper)
#> Y1 Y2 Y3 Y1 Y2 Y3 Y1
#> Y1 0.04284 0.00544 -0.01684 . 0.03780 -0.01297 -0.04097 . 0.04842
#> Y2 0.00544 1.31266 0.06543 . -0.01297 1.15443 -0.06808 . 0.02483
#> Y3 -0.01684 0.06543 1.84833 . -0.04097 -0.06808 1.62191 . 0.00501
#> Y2 Y3
#> Y1 0.02483 0.00501
#> Y2 1.48829 0.18637
#> Y3 0.18637 2.07349
#>
#>
#> Regime2:
#>
#> Autoregressive coefficients
#> Mean 2(1-PD) HDI.Lower HDI.Upper Mean 2(1-PD)
#> (Intercept) 0.19806 1e-05 0.17236 0.22438 | 0.19580 1e-05
#> Y1.lag(1) 0.22521 1e-05 0.19085 0.26044 | -0.29698 1e-05
#> Y2.lag(1) -0.11839 1e-05 -0.14077 -0.09683 | 0.19611 1e-05
#> Y3.lag(1) -0.08728 1e-05 -0.09856 -0.07596 | 0.32308 1e-05
#> Y1.lag(2) -0.13281 1e-05 -0.16558 -0.09741 | -0.16747 1e-05
#> Y2.lag(2) -0.34440 1e-05 -0.36325 -0.32421 | 0.19027 1e-05
#> Y3.lag(2) 0.21655 1e-05 0.20225 0.23019 | 0.16870 1e-05
#> HDI.Lower HDI.Upper Mean 2(1-PD) HDI.Lower HDI.Upper
#> (Intercept) 0.15887 0.23039 | -0.24841 0.00001 -0.34754 -0.14078
#> Y1.lag(1) -0.34693 -0.24798 | 0.26756 0.00067 0.12847 0.41148
#> Y2.lag(1) 0.16790 0.22617 | 0.13642 0.00200 0.04436 0.21669
#> Y3.lag(1) 0.30738 0.33920 | -0.23078 0.00001 -0.27636 -0.18479
#> Y1.lag(2) -0.21448 -0.11972 | -0.29378 0.00001 -0.42693 -0.15278
#> Y2.lag(2) 0.16187 0.21604 | -0.14354 0.00067 -0.22095 -0.06226
#> Y3.lag(2) 0.14956 0.18885 | 0.21172 0.00001 0.15407 0.26757
#>
#> Scale parameter (Mean, HDI.Lower, HDI.Upper)
#> Y1 Y2 Y3 Y1 Y2 Y3 Y1
#> Y1 0.14711 -0.00924 0.01547 . 0.13275 -0.02183 -0.02059 . 0.16166
#> Y2 -0.00924 0.29331 -0.04912 . -0.02183 0.26423 -0.10076 . 0.00320
#> Y3 0.01547 -0.04912 2.46329 . -0.02059 -0.10076 2.23850 . 0.05194
#> Y2 Y3
#> Y1 0.00320 0.05194
#> Y2 0.32189 0.00017
#> Y3 0.00017 2.70172
#>
#>
#> Extra parameter
#> Mean 2(1-PD) HDI.Lower HDI.Upper
#> nu 6.24436 . 5.4349 7.21871
#>
#>
###### Simulation of a trivariate VAR model
n <- 2000
k <- 3
myars <- ars(nregim=1,p=2)
dist <- "Slash"; extra <- 2
parms <- list()
for(j in 1:myars$nregim){
np <- 1 + myars$p[j]*k
parms[[j]] <- list()
parms[[j]]$location <- c(ifelse(runif(np*k)<=0.5,1,-1)*rbeta(np*k,shape1=4,shape2=16))
parms[[j]]$location <- matrix(parms[[j]]$location,np,k)
parms[[j]]$scale <- rgamma(k,shape=1,scale=1)*diag(k)
}
out2 <- simtar(n=n, k=k, ars=myars, parms=parms, dist=dist, extra=extra)
#>
#>
#> Sample size :2000 time points
#>
#> Output Series :Y1 | Y2 | Y3
#>
#> Error Distribution :Slash(2)
#>
#> Number of regimes :1
#>
#> Deterministics :Intercept
#>
#> Autoregressive orders:2 in each regime
#>
#>
#> Regime 1:
#>
#> Autoregressive coefficients
#> Intercept phi_1 phi_2
#> Y1 -0.1328379 | -0.2120508 0.1259517 -0.1305420 | -0.1166168 0.1202440
#> Y2 -0.1494702 | -0.1659424 -0.3989037 0.3255796 | 0.1188607 0.3324750
#> Y3 0.2215310 | 0.2388621 -0.1116290 0.1471736 | -0.1133576 0.1261956
#>
#> Y1 0.1079300 |
#> Y2 -0.2017739 |
#> Y3 0.3369314 |
#>
#> Scale parameter
#> Y1 Y2 Y3
#> Y1 0.7463516 0.000000 0.0000000
#> Y2 0.0000000 1.666729 0.0000000
#> Y3 0.0000000 0.000000 0.3769205
str(out2)
#> 'data.frame': 2002 obs. of 4 variables:
#> $ Y1 : num 1.0134 1.0627 -0.0337 -0.8272 3.609 ...
#> $ Y2 : num -0.539 0.679 3.272 1.682 -1.244 ...
#> $ Y3 : num 0.921 0.583 3.497 0.734 1.781 ...
#> $ Regime: num NA NA 1 1 1 1 1 1 1 1 ...
fit2 <- mtar(~ Y1 + Y2 + Y3, data=out2, ars=myars, dist=dist,
n.burn=2000, n.sim=3000, n.thin=2)
summary(fit2)
#>
#>
#> Sample size : 2000 time points
#>
#> Output Series : Y1 | Y2 | Y3
#>
#> Error Distribution : Slash
#>
#> Number of regimes : 1
#>
#> Deterministics : Intercept
#>
#> Autoregressive orders: 2 in each regime
#>
#>
#>
#>
#> Regime1:
#>
#> Autoregressive coefficients
#> Mean 2(1-PD) HDI.Lower HDI.Upper Mean 2(1-PD)
#> (Intercept) -0.14877 1e-05 -0.21144 -0.09160 | -0.05178 0.26400
#> Y1.lag(1) -0.21400 1e-05 -0.23190 -0.19755 | -0.15526 0.00001
#> Y2.lag(1) 0.12629 1e-05 0.11299 0.14091 | -0.39226 0.00001
#> Y3.lag(1) -0.15383 1e-05 -0.18319 -0.12441 | 0.30089 0.00001
#> Y1.lag(2) -0.10207 1e-05 -0.11681 -0.08577 | 0.13268 0.00001
#> Y2.lag(2) 0.11334 1e-05 0.09754 0.12833 | 0.32896 0.00001
#> Y3.lag(2) 0.11624 1e-05 0.09268 0.14045 | -0.17817 0.00001
#> HDI.Lower HDI.Upper Mean 2(1-PD) HDI.Lower HDI.Upper
#> (Intercept) -0.15027 0.04013 | 0.20854 1e-05 0.16334 0.25268
#> Y1.lag(1) -0.17935 -0.13253 | 0.23496 1e-05 0.22050 0.24714
#> Y2.lag(1) -0.41323 -0.37121 | -0.10841 1e-05 -0.11821 -0.09885
#> Y3.lag(1) 0.25590 0.34759 | 0.15607 1e-05 0.13518 0.17745
#> Y1.lag(2) 0.10780 0.15703 | -0.11472 1e-05 -0.12695 -0.10398
#> Y2.lag(2) 0.30543 0.35149 | 0.12782 1e-05 0.11642 0.13863
#> Y3.lag(2) -0.21431 -0.14037 | 0.33523 1e-05 0.31724 0.35169
#>
#> Scale parameter (Mean, HDI.Lower, HDI.Upper)
#> Y1 Y2 Y3 Y1 Y2 Y3 Y1
#> Y1 0.69594 0.00473 0.00786 . 0.63374 -0.05024 -0.01659 . 0.76310
#> Y2 0.00473 1.69761 -0.03083 . -0.05024 1.53347 -0.07172 . 0.05510
#> Y3 0.00786 -0.03083 0.37444 . -0.01659 -0.07172 0.34006 . 0.03335
#> Y2 Y3
#> Y1 0.05510 0.03335
#> Y2 1.85559 0.00674
#> Y3 0.00674 0.41039
#>
#>
#> Extra parameter
#> Mean 2(1-PD) HDI.Lower HDI.Upper
#> nu 2.0093 . 1.86607 2.15787
#>
#>
n <- 5000
k <- 3
myars <- ars(nregim=2,p=c(1,2))
dist <- "Laplace"
parms <- list()
for(j in 1:myars$nregim){
np <- 1 + myars$p[j]*k
parms[[j]] <- list()
parms[[j]]$location <- c(ifelse(runif(np*k)<=0.5,1,-1)*rbeta(np*k,shape1=4,shape2=16))
parms[[j]]$location <- matrix(parms[[j]]$location,np,k)
parms[[j]]$scale <- rgamma(k,shape=1,scale=1)*diag(k)
}
out3 <- simtar(n=n, k=k, ars=myars, parms=parms, delay=2,
thresholds=-1, dist=dist, setar=2)
#>
#>
#> Sample size :5000 time points
#>
#> Output Series :Y1 | Y2 | Y3
#>
#> Threshold Series :Y2 with a delay equal to 2
#>
#> Error Distribution :Laplace
#>
#> Number of regimes :2
#>
#> Deterministics :Intercept
#>
#> Autoregressive orders:1, 2
#>
#>
#> Thresholds
#>
#> Regime 1 (-Inf,-1]
#> Regime 2 (-1,Inf)
#>
#>
#> Regime 1:
#>
#> Autoregressive coefficients
#> Intercept phi_1
#> Y1 -0.3955847 | -0.25919156 -0.08655823 -0.2749034 |
#> Y2 -0.2756847 | 0.09024662 0.21837686 -0.1522644 |
#> Y3 -0.1182234 | 0.12287218 0.07297322 0.2214600 |
#>
#> Scale parameter
#> Y1 Y2 Y3
#> Y1 1.602306 0.000000 0.000000
#> Y2 0.000000 2.162429 0.000000
#> Y3 0.000000 0.000000 3.244577
#>
#> Regime 2:
#>
#> Autoregressive coefficients
#> Intercept phi_1 phi_2
#> Y1 -0.30734113 | 0.2223958 0.12902632 0.2476028 | -0.1839958 0.1430500
#> Y2 -0.09297533 | -0.1449571 -0.06225941 -0.3168271 | -0.1277287 -0.2456443
#> Y3 -0.26178469 | -0.2618776 0.17515898 -0.1218319 | -0.1445789 0.2048370
#>
#> Y1 -0.1676492 |
#> Y2 0.1221158 |
#> Y3 -0.2190242 |
#>
#> Scale parameter
#> Y1 Y2 Y3
#> Y1 0.2349453 0.0000000 0.00000000
#> Y2 0.0000000 0.1437912 0.00000000
#> Y3 0.0000000 0.0000000 0.04735645
str(out3)
#> 'data.frame': 5002 obs. of 4 variables:
#> $ Y1 : num 1.309 -0.687 -0.618 0.443 -0.725 ...
#> $ Y2 : num 1.928 0.93 0.116 -2.075 -0.432 ...
#> $ Y3 : num -0.902 -1.058 0.243 0.779 -0.4 ...
#> $ Regime: num NA NA 2 2 2 1 2 2 2 2 ...
fit3 <- mtar(~ Y1 + Y2 + Y3, data=out3, ars=myars, dist=dist,
n.burn=2000, n.sim=3000, n.thin=2, setar=2)
summary(fit3)
#>
#>
#> Sample size : 4999 time points
#>
#> Output Series : Y1 | Y2 | Y3
#>
#> Threshold Series : Y2 with a estimated delay equal to 2
#>
#> Error Distribution : Laplace
#>
#> Number of regimes : 2
#>
#> Deterministics : Intercept
#>
#> Autoregressive orders: 1, 2
#>
#>
#>
#> Thresholds (Mean, HDI.Lower, HDI.Upper)
#>
#> Regime 1 (-Inf,-0.29515] (-Inf,-0.3069] (-Inf,-0.25855]
#> Regime 2 (-0.29515,Inf) (-0.3069,Inf) (-0.25855,Inf)
#>
#>
#> Regime1:
#>
#> Autoregressive coefficients
#> Mean 2(1-PD) HDI.Lower HDI.Upper Mean 2(1-PD)
#> (Intercept) -0.37433 0.00001 -0.44444 -0.31104 | -0.03319 0.31600
#> Y1.lag(1) -0.01744 0.44467 -0.06515 0.02575 | -0.02073 0.32667
#> Y2.lag(1) -0.01072 0.58467 -0.05006 0.02689 | 0.04337 0.02800
#> Y3.lag(1) 0.00435 0.86800 -0.04501 0.04832 | -0.19588 0.00001
#> HDI.Lower HDI.Upper Mean 2(1-PD) HDI.Lower HDI.Upper
#> (Intercept) -0.09671 0.03091 | -0.29534 0.00001 -0.36027 -0.22468
#> Y1.lag(1) -0.06190 0.02011 | -0.11110 0.00001 -0.15588 -0.06526
#> Y2.lag(1) 0.00579 0.08202 | 0.16773 0.00001 0.13243 0.20276
#> Y3.lag(1) -0.23031 -0.16075 | -0.01024 0.60533 -0.04933 0.02922
#>
#> Scale parameter (Mean, HDI.Lower, HDI.Upper)
#> Y1 Y2 Y3 Y1 Y2 Y3 Y1
#> Y1 1.30380 -0.11629 -0.07522 . 1.20845 -0.17671 -0.14852 . 1.41401
#> Y2 -0.11629 1.31929 0.05850 . -0.17671 1.21325 -0.00655 . -0.05797
#> Y3 -0.07522 0.05850 1.84404 . -0.14852 -0.00655 1.69281 . -0.00598
#> Y2 Y3
#> Y1 -0.05797 -0.00598
#> Y2 1.41967 0.13334
#> Y3 0.13334 1.99726
#>
#>
#> Regime2:
#>
#> Autoregressive coefficients
#> Mean 2(1-PD) HDI.Lower HDI.Upper Mean 2(1-PD)
#> (Intercept) -0.28137 1e-05 -0.31927 -0.24552 | -0.11014 1e-05
#> Y1.lag(1) 0.22868 1e-05 0.21573 0.24161 | -0.14412 1e-05
#> Y2.lag(1) 0.12908 1e-05 0.11700 0.14050 | -0.06099 1e-05
#> Y3.lag(1) 0.25299 1e-05 0.24191 0.26271 | -0.31581 1e-05
#> Y1.lag(2) -0.18214 1e-05 -0.19178 -0.17112 | -0.12984 1e-05
#> Y2.lag(2) 0.13356 1e-05 0.12160 0.14610 | -0.24129 1e-05
#> Y3.lag(2) -0.16933 1e-05 -0.17828 -0.15999 | 0.11627 1e-05
#> HDI.Lower HDI.Upper Mean 2(1-PD) HDI.Lower HDI.Upper
#> (Intercept) -0.13517 -0.08506 | -0.25176 1e-05 -0.26768 -0.23755
#> Y1.lag(1) -0.15377 -0.13405 | -0.25566 1e-05 -0.26099 -0.24939
#> Y2.lag(1) -0.07171 -0.05092 | 0.17967 1e-05 0.17461 0.18558
#> Y3.lag(1) -0.32401 -0.30718 | -0.12243 1e-05 -0.12713 -0.11796
#> Y1.lag(2) -0.13693 -0.12076 | -0.14533 1e-05 -0.14944 -0.14087
#> Y2.lag(2) -0.25168 -0.23135 | 0.20066 1e-05 0.19476 0.20604
#> Y3.lag(2) 0.10758 0.12552 | -0.21709 1e-05 -0.22153 -0.21308
#>
#> Scale parameter (Mean, HDI.Lower, HDI.Upper)
#> Y1 Y2 Y3 Y1 Y2 Y3 Y1
#> Y1 0.22043 -0.00335 -0.00041 . 0.20344 -0.01094 -0.00487 . 0.23702
#> Y2 -0.00335 0.13979 0.00037 . -0.01094 0.12970 -0.00336 . 0.00487
#> Y3 -0.00041 0.00037 0.04580 . -0.00487 -0.00336 0.04240 . 0.00411
#> Y2 Y3
#> Y1 0.00487 0.00411
#> Y2 0.15110 0.00387
#> Y3 0.00387 0.04931
#>
#>
# }