EViews.com

Hi all,

I have three questions on STAR output.

I've estimated a univariate ESTAR model with 4 lags on both linear and non linear parts and a delay of 2 in both EVIEWS and R Studio. All coefficients are the same apart from gamma (R Studio gamma is 1.1359, std error 0.3610, sigma2 of model 23.256613). I know there is an issue with standardising gamma but I think it should be possile to convert by hand by scaling with the std dev fof the dep var, but nothing seems to come close? Can anyone help here?

When performing the linearity tests some of the p values are reported as NA. Can anyone explain this?
Also, with the same tests, the f statistics are the same for several hypothesis tests. Why is this?

All relevant details in the attachments below.

Thanks.

You'll need to post your workfile.

Thanks, here it is.

Which R package are you using?

RStudio R version 4.3.0 (2023-04-21 ucrt) -- "Already Tomorrow"

Package RSTAR.

Actually if you know another package that does exponential STAR models I'd love to hear it too. There's a few that does the LSTAR but as far as I know RSTAR is the only ESTAR package. It's great but there's no help file available.

Hi Rebecca,

Also, in trying to estimate a ESTAR with lagged differences -something like [y @base(c) @base(d(y)(-1) @alt(y(-1)) @alt(d(y)(-1))] or any other effort to introduce lagged differences (a common version of the STAR model) EVIEWS crashes without error message.

Not sure how to get around this and the ohter three issue mentioned above. Any help greatly appreciated.

Steve

I was asking because I don't know of any R package that does ESTAR.

By RSTAR, do you mean https://github.com/mbalcilar/RSTAR? Unfortunately, this doesn't appear to include any code. Without access to the underlying R code I can't see how they're doing their gamma calculation. On your end, however, you can try one of the many methods for examining the insides of R packages, such as typing the name of the function at the prompt.

Hi Rebecca,

Also, in trying to estimate a ESTAR with lagged differences -something like [y @base(c) @base(d(y)(-1) @alt(y(-1)) @alt(d(y)(-1))] or any other effort to introduce lagged differences (a common version of the STAR model) EVIEWS crashes without error message.

Not sure how to get around this and the ohter three issue mentioned above. Any help greatly appreciated.

Steve

Are you trying to do d(y(-1)), etc?

Thanks Rebecca, the brackets were causing it to crash. Stupid of me.

Can you comment on the NA issue and the F stats?

Here is the function code. Understanding it is beyond me I'm afraid.

function (x, p1, p2, d = NULL, linear = NULL, transition = NULL,
intercept = TRUE, model = c("LSTAR", "ESTAR"), g.scale, hccvcov = FALSE,
title = NULL, description = NULL, lab.transition = NULL,
lab.x = NULL, param0 = NULL, grid.start = TRUE, a.grad = TRUE,
gamma.upper = 200, gamma.lower = exp(-4.99), single.p1 = FALSE,
single.p2 = FALSE, ...)
{
model = match.arg(model)
if (missing(p1) && missing(p1))
stop("missing p1 and p2")
if (is.null(transition) && is.null(d))
stop("specify transition or d")
if (!is.null(transition) && !is.null(d))
stop("specify only one of transition and d")
if (p1 == 0 && p2 == 0)
stop("both p1 and p2 are 0")
if (!is.null(linear) && intercept)
linear = cbind(1, linear)
if (is.null(linear)) {
p = 0
k0 = 0
}
else {
linear = as.matrix(linear)
k0 = ncol(linear)
}
if (missing(p1)) {
p1 = 0
k1 = 0
}
else if (length(p1) == 1) {
if (p1 != 0) {
if (single.p1) {
k1 = 1
}
else {
p1 = 1:p1
k1 = length(p1)
}
}
else {
k1 = 0
}
}
else k1 = length(p1)
if (missing(p2)) {
p2 = 0
k2 = 0
}
else if (length(p2) == 1) {
if (p2 != 0) {
if (single.p2) {
k2 = 1
}
else {
p2 = 1:p2
k2 = length(p2)
}
}
else {
k2 = 0
}
}
else k2 = length(p2)
if (k0 > 0) {
sel0 = 1:k0
sk0 = length(sel0)
p.add = 0
}
else {
sel0 = 0
sk0 = 0
p.add = 1
}
if (k1 > 0) {
sel1 = (sk0 + 1):(sk0 + k1 + p.add)
sk1 = length(sel1)
}
else {
sel1 = 0
sk1 = 0
}
if (k2 > 0) {
sel2 = (sk0 + sk1 + 1):(sk0 + sk1 + k2 + p.add)
sk2 = length(sel2)
}
else {
sel2 = 0
sk2 = 0
}
par.sel = list(sel0 = sel0, sel1 = sel1, sel2 = sel2, sk0 = sk0,
sk1 = sk1, sk2 = sk2)
if (!is.null(transition)) {
if (!is.null(d))
warning("ignoring d")
d = NULL
z = transition
z.name = deparse(substitute(transition))
}
x.name = deparse(substitute(x))
p.max = max(c(max(p1), max(p2)))
if (!is.null(d)) {
if (p.max > d)
data <- embed(x, p.max + 1)
else data <- embed(x, d + 1)
}
else data <- embed(x, p.max + 1)
k = k0 + k1 + k2
xo = x
nxo = length(xo)
ists <- is.ts(x)
if (!ists)
x <- as.ts(x)
x.freq <- frequency(x)
if (any(is.na(x)))
stop("NAs in x")
x.tsp <- tsp(x)
x = as.matrix(data[, -1])
y = data[, 1]
ny = NROW(x)
nx = nrow(x)
if (k1 > 0) {
if (k0 > 0)
x1 = x[, p1]
else x1 = cbind(1, x[, p1])
}
else x1 = numeric()
if (k2 > 0) {
if (k0 > 0)
x2 = x[, p2]
else x2 = cbind(1, x[, p2])
}
else x2 = numeric()
x = cbind(x1, x2)
if (!is.null(d)) {
z = data[, d + 1]
z.name = paste(x.name, "(t-", d, ")", sep = "")
}
else {
nz <- NROW(z)
adj = nz - ny
z = c(z)[(adj + 1):nz]
}
if (k0 > 0) {
nl = nrow(linear)
adj = nl - nx
x = cbind(linear[(adj + 1):nl, ], x)
}
n = NROW(y)
if (missing(param0))
param0 = NULL
if (missing(g.scale)) {
if (model == "LSTAR")
g.scale = sd(z)
else g.scale = var(z)
}
if (model == "ESTAR")
fit = ESTAR(x, y, z, par.sel, g.scale, grid.start, a.grad,
param0, hccvcov, gamma.upper, gamma.lower)
else if (model == "LSTAR")
fit = LSTAR(x, y, z, par.sel, g.scale, grid.start, a.grad,
param0, hccvcov, gamma.upper, gamma.lower)
else stop("uknown model type")
G = fit$G
par = fit$coef
knl <- length(par) - (sk0 + sk1 + sk2)
X = getxmat.lin(x, par.sel, G)
if (par.sel$sk0 > 0) {
fit.linear = c(X[, par.sel$sel0, drop = FALSE] %*% par[(knl +
1):(knl + sk0)])
fit.linear = ts(fit.linear, frequency = x.tsp[3], end = x.tsp[2])
}
else {
fit.linear = NULL
}
if (any(p1 > 0)) {
fit.trs1 = c(X[, par.sel$sel1, drop = FALSE] %*% par[(knl +
sk0 + 1):(knl + sk0 + sk1)])
fit.trs1 = ts(fit.trs1, frequency = x.tsp[3], end = x.tsp[2])
}
else {
fit.trs1 = NULL
}
if (any(p2 > 0)) {
fit.trs2 = c(X[, par.sel$sel2, drop = FALSE] %*% par[(knl +
sk0 + sk1 + 1):(knl + sk0 + sk1 + sk2)])
fit.trs2 = ts(fit.trs2, frequency = x.tsp[3], end = x.tsp[2])
}
else {
fit.trs2 = NULL
}
if (model == "LSTAR")
par.names = c("Gamma", "Treshold")
if (model == "LSTR2")
par.names = c("Gamma", "Treshold_1", "Treshold_2")
if (model == "ESTAR")
par.names = c("Gamma", "Treshold_1")
if (model == "ESTR2")
par.names = c("Gamma", "Treshold_1", "Treshold_2")
if (k0 > 0)
if (intercept)
par.names = c(par.names, "Intercep", paste("X0_",
1:(k0 - 1), sep = ""))
else par.names = c(par.names, paste("X0_", 1:k0, sep = ""))
if (k1 > 0) {
if (k0 > 0)
par.names = c(par.names, paste("AR_1(", p1, ")",
sep = ""))
else par.names = c(par.names, "Intercep_1", paste("AR_1(",
p1, ")", sep = ""))
}
if (k2 > 0) {
if (k0 > 0)
par.names = c(par.names, paste("AR_2(", p2, ")",
sep = ""))
else par.names = c(par.names, "Intercep_2", paste("AR_2(",
p2, ")", sep = ""))
}
k = length(fit$coef)
s2 = fit$sigma2
logL = -(n/2) * log(2 * pi) - (n/2) * log(s2) - (1/(2 * s2)) *
(n * s2)
inf.crit = info.crit(logL, k, n)
fit$info.crit = inf.crit
names(fit$coef) = par.names
colnames(fit$gradient) = par.names
fit$se.coef = fit$se
fit$var.coef
fit$tstitle = paste(model, "(p1,p2) with transition variable = ",
z.name, sep = "")
fit$tsmodel = model
fit$linear = fit.linear
fit$transition1 = fit.trs1
fit$transition2 = fit.trs2
fit$class = "STAR"
class(fit) = "list"
if (is.null(title))
title = "STAR Modelling"
if (is.null(description))
description = as.character(date())
if (is.null(lab.transition))
lab.transition = z.name
if (is.null(lab.x))
lab.x = x.name
residuals = c(fit$residuals)
residuals = ts(residuals, frequency = x.tsp[3], end = x.tsp[2])
fitted = c(fit$fitted.values)
fitted = ts(fitted, frequency = x.tsp[3], end = x.tsp[2])
y = ts(y, frequency = x.tsp[3], end = x.tsp[2])
fit$residuals <- residuals
fit$fitted.values <- fitted
new("STAR", call = as.call(match.call()), model = as.character(model),
parameter = list(p1 = p1, p2 = p2, d = d, transition = z,
lab.x = lab.x, lab.transition = lab.transition, par.sel = par.sel,
knl = knl, g.scale = g.scale, n.used = n, grid.start = grid.start,
a.grad = a.grad, gamma.upper = gamma.upper, gamma.lower = gamma.lower,
param0 = param0), data = list(x = y, x.reg = x, x.org = xo),
fit = fit, residuals = residuals, fitted.values = fitted,
title = as.character(title), description = as.character(description))
}
<bytecode: 0x000001faa6fd35a8>
<environment: namespace:RSTAR>

Partway down it calls the ESTAR function. Look inside.

Ok. I’m away from the machine for a week now but I’ll post on Tuesday week. Thanks.

For the F stats: we're looking into it. Current thinking is there probably aren't enough observations for some of the tests.

Hi Rebecca

> ESTAR
function (x, y, z, par.sel, g.scale = var(z), grid.start = TRUE,
a.grad = TRUE, param0 = NULL, hccvcov = FALSE, gamma.upper = 200,
gamma.lower = exp(-4.99))
{
if (gamma.lower < 1e-04) {
gamma.lower <- exp(-4.99)
warning("lower limit of gamma too small, setting it to exp(-4.99)")
}
if (grid.start)
start = startstar(x = x, y = y, z = z, par.sel = par.sel,
g.scale = g.scale, model = "ESTAR", param0 = param0,
gamma.upper = gamma.upper, gamma.lower = gamma.lower)
else {
if (is.null(param0)) {
start = c(0, mean(z))
}
else start = param0
}
n = length(y)
s.z = sort(z)
lower = c(log(gamma.lower), s.z[ceiling(0.05 * n)])
upper = c(log(gamma.upper), s.z[floor(0.95 * n)])
grd <- function(gamc, ...) {
res <- ESTAR.res(gamc, ...)
theta <- c(gamc, res$ar)
dgamdc <- dexponential(theta, x, y, z, par.sel, g.scale)
grad = -2 * (t(dgamdc) %*% res$residual)
return(grad)
}
if (a.grad) {
o <- optim(par = start, fn = ESTAR.ssr, gr = grd, lower = lower,
upper = upper, method = "L-BFGS-B", hessian = TRUE,
control = list(), y = y, x = x, z = z, par.sel = par.sel,
g.scale = g.scale)
}
else {
o <- optim(par = start, fn = ESTAR.ssr, gr = NULL, lower = lower,
upper = upper, method = "L-BFGS-B", hessian = TRUE,
control = list(), y = y, x = x, z = z, par.sel = par.sel,
g.scale = g.scale)
}
res <- ESTAR.res(o$par, x, y, z, par.sel, g.scale)
par = c(o$par, ar = res$ar)
varcov = vcovestar(par, x, y, z, c(res$residual), par.sel,
g.scale)
if (hccvcov)
se2 = sqrt(diag(varcov$hccvcov))
else se2 = sqrt(diag(varcov$vcov))
gradient <- gradestar(par, x, y, z, par.sel, g.scale)
n = length(y)
residual <- res$residual
fitted <- y - res$residual
G = exponential(z, g.scale, exp(par[1]), par[2])
par[1] = exp(par[1])
sigma2 = o$value/n
return(list(coef = par, se.coef = se2, sigma2 = sigma2, residuals = residual,
fitted.values = fitted, G = G, var.coef = varcov$vcov,
var.coef.hcc = varcov$hccvcov, gradient = gradient, convergence = o$convergence,
message = o$message))
}
<bytecode: 0x0000017f34e1b4c8>
<environment: namespace:RSTAR>

There's a lot going on in these code snippets. As a start, can you use trace in R (trace("functionname", edit=TRUE)) to test how changing the log transformation of gamma in ESTAR affects your results?

Ok thanks Rebecca. It will be a week or so before I can post this.