1 |
z.ts |
|
etw |
|
p |
|
q |
|
n |
|
r |
|
type |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 | ##---- Should be DIRECTLY executable !! ----
##-- ==> Define data, use random,
##-- or do help(data=index) for the standard data sets.
## The function is currently defined as
function(z.ts,etw,p,q=NULL,n=NULL,r,type="vecm")
# z.ts ... time series
# etw ... estimation time window as specified in estim.we.mdls.r
# p ... lag order
# n ... number of endogenous variables
# r ... cointegrating relations
# type ... model, either "vecm" or "we.vecm"
{
if (is.null(n) && type=="we.vecm") stop("\'n\' missing",call.=FALSE)
mdls <- list()
test.stat <- vector()
cv <- vector()
p.value <- vector()
T <- dim(z.ts)[1]-p
m <- dim(z.ts)[2]
cases <- c("I","II","III","IV","V")
if (type=="we.vecm")
{
# Johansen ~ Pesaran
# case<- "H_2(r)" ~ case "I" # \mu_t = 0
# case<- "H_1^*(r)" ~ case "II" # \mu_t = \alpha \rho_0
# case<- "H_1(r)" ~ case "III" # \mu_t = \mu_0
# case<- "H^*(r)" ~ case "IV" # \mu_t = \mu_0+\alpha \rho_1 t # \alpha_{\bot}'\mu_1=0
# case<- "H(r)" ~ case "V" # \mu_t = \mu_0+\mu_1 t
for (i in 1:5){mdls[[i]] <- est.we.mdls(z.ts=z.ts,etw=etw,p=p,q=q,r=r,n=n,case=cases[i])}
} else {
for (i in 1:5){mdls[[i]] <- est.vecm.mdls(Y.ts=z.ts,etw=etw,p=p,r=r,case=cases[i])}
}
# test statistics as in Johansen (1995) p. 161 f.
for (i in c(1,3))
{
test.stat[i] <- T*sum(log(1-mdls[[5-i+1]]$lambda[(r+1):m])-log(1-mdls[[4-i+1]]$lambda[(r+1):m]))
cv[i] <- qchisq(0.95,df=(m-r))
p.value[i] <- 1-pchisq(test.stat[i],df=(m-r))
}
for (i in c(2,4))
{
test.stat[i] <- T*sum(log(1-mdls[[3-i+2]]$lambda[1:r])-log(1-mdls[[4-i+2]]$lambda[1:r]))
cv[i] <- qchisq(0.95,df=r)
p.value[i] <- 1-pchisq(test.stat[i],df=r)
}
# prepare & print output
res <- cbind(round(test.stat,4),round(cv,4),round(p.value,4))
rownames(res) <- c("c.1 restricted vs. c.1 unrestricted"," c.1 = 0 vs. c.1 restricted",
"c.0 restricted vs. c.0 unrestricted"," c.0 = 0 vs. c.0 restricted")
colnames(res) <- c(" Statistic"," Crit.value"," p-value")
cat("\nCase selection results:\n")
print(res)
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.