您好！我太菜了，问题不断！如果想估计下面（图片格式）这个模型，其中G(.)为logistic函数，在R中该如何编程序呢？谢谢您了啊！
30# epoh

LSTAR model基本结构不变

G <- function(y, g, th)
1/(1 + exp(-g*(y-th)))


F <- function(phi1, phi2, g, th){

xxL %*% phi1 + (xxH %*% phi2) * G(z, g, th)}
只是数据有所改变
以你新的new model,log10(lynx),m=2,d=1 为例
       y(t-1)      dy(t-1)       dy(t)
  [1,] 2.506505  0.076752752  0.260650834
  [2,] 2.767156  0.260650834  0.172862289
  [3,] 2.940018  0.172862289  0.228773865
       ...........
log10(lynx),m=3,d=1
        y(t-1)      dy(t-2)     dy(t-1)       dy(t)
  [1,] 2.767156  0.076752752  0.260650834  0.172862289
  [2,] 2.940018  0.260650834  0.172862289  0.228773865
  [3,] 3.168792  0.172862289  0.228773865  0.281611066

lstar_newmodel.R

lstar_newmodel.rar
大小:(6.16 KB)

 马上下载

本附件包括：

• lstar_newmodel.R

source("lstar_newmodel.R")
mod.lstar <- lstar(log10(lynx), m=2, d=1, control=list(maxit=3000))
mod.lstar
summary(mod.lstar)

mod1.lstar <- lstar(log10(lynx), m=3, d=1, control=list(maxit=3000))
mod1.lstar
summary(mod1.lstar)

LSTAR model
Coefficients:
Low regime:
    phi1.0     phi1.1     phi1.2     phi1.3
0.5389244 -0.1335123  0.1276239  0.2615295
High regime:
     phi2.0      phi2.1      phi2.2      phi2.3
0.83454679 -0.34393225  0.09368726  0.63782175
Smoothing parameter: gamma = 40.03
Threshold
Variable: Z(t) = + (1) X(t) + (0) X(t-1)+ (0) X(t-2)
Value: 2.568
Residuals:
      Min        1Q    Median        3Q       Max
-0.577726 -0.116138  0.025577  0.140332  0.439914
Fit:
residuals variance = 0.03941,  AIC = -349, MAPE = 128.3%
Non-linearity test of full-order LSTAR model against full-order AR model
F = 5.8493 ; p-value = 0.00098792
Threshold
Variable: Z(t) = + (1) X(t) + (0) X(t-1)+ (0) X(t-2)

#######cpi
data=read.table("data.txt", header = TRUE)
x <- data$cpi
y <- data$rpi
z <- data$m2
cpi.lstar <- lstar(x, m=3,d=1, thVar=z,control=list(maxit=3000))
cpi.lstar
summary(cpi.lstar)
LSTAR model
Coefficients:
Low regime:
     phi1.0      phi1.1      phi1.2      phi1.3
-0.28408868 -0.04046736 -0.49254774  0.97649199
High regime:
   phi2.0    phi2.1    phi2.2    phi2.3
4.397429 -1.083845  0.993229 -1.255000
Smoothing parameter: gamma = 48.02
Threshold
Variable: external
Value: 0.3224
Residuals:
     Min       1Q   Median       3Q      Max
-5.52572 -1.58668 -0.30896  1.65844  4.79960
Fit:
residuals variance = 6.057,  AIC = 54, MAPE = 168.2%
Non-linearity test of full-order LSTAR model against full-order AR model
F = 0.11978 ; p-value = 0.94616
Threshold
Variable: external

您好！再向您请教个问题，我用cr=cpi.lstar$residuals调用了残差后，如果我想对残差进行独立性检验（Q统计量）、正态性检验（JB统计量）和ARCH检验（LM统计量，也可以是其他统计量），该如何编程呢？谢谢您了！ 32# epoh

前辈您好，转换变量是滞后变量时我可以用S-PLUS做出转换变量和转换函数的散点图，但是如果转换变量是外生变量时我只能用R估计模型，但是该如何调用转换函数的值呢？我试了 rff=rdx.lstar$transition.func> rff
NULL
可是不能调用，希望您能解答，谢谢您了！ 21# epoh

epoch是很好的老师，非常感谢

source("lstar_newmodel.R")
data=read.table("data.txt", header = TRUE)
x <- data$cpi
y <- data$rpi
z <- data$m2
cpi.lstar <- lstar(x, m=2,d=1, thVar=z,control=list(maxit=3000))
cpi.lstar
names(cpi.lstar)
#[1] "str"      "coefficients"   "fitted.values"  "residuals"   "k"            
#[6] "model"    "model.specific"
g=cpi.lstar$coefficients[7]
th=cpi.lstar$coefficients[8]
gfun <- function(y, g, th)      1 / (1 + exp(-g*(y-th)))  
G=gfun(z,g,th)
par(mfrow=c(2,1))
plot(G,type = "l", col = "red")
plot(G,z)

#########Statistical Tests:
library(fGarch)
# Lagged Series:
    .tslagGarch = function (x, k = 1) {
        ans = NULL
        for (i in k) ans = cbind(ans, .tslag1Garch(x, i))
        indexes = (1:length(ans[, 1]))[!is.na(apply(ans, 1, sum))]
        ans = ans[indexes, ]
        if (length(k) == 1) ans = as.vector(ans)
        ans }
    .tslag1Garch = function (x, k) {
        c(rep(NA, times = k), x[1:(length(x) - k)]) }
    # Statistical Tests:
    cat("\n Residuals Tests:\n")
    r.s = cpi.lstar$residuals
    ans = NULL
    # Normality Tests:
    jbtest = jarqueberaTest(r.s)@test
    ans = rbind(ans, c(jbtest[1], jbtest[2]))
    if (length(r.s) < 5000) {
        swtest = shapiro.test(r.s)
        if (swtest[2] < 2.6e-16) swtest[2] = 0
        ans = rbind(ans, c(swtest[1], swtest[2]))
    } else {
        ans = rbind(ans, c(NA, NA))
    }
    # Ljung-Box Tests:
    box10 = Box.test(r.s, lag = 10, type = "Ljung-Box")
    box15 = Box.test(r.s, lag = 15, type = "Ljung-Box")
    box20 = Box.test(r.s, lag = 20, type = "Ljung-Box")
    ans = rbind(ans, c(box10[1], box10[3]))
    ans = rbind(ans, c(box15[1], box15[3]))
    ans = rbind(ans, c(box20[1], box20[3]))
    box10 = Box.test(r.s^2, lag = 10, type = "Ljung-Box")
    box15 = Box.test(r.s^2, lag = 15, type = "Ljung-Box")
    box20 = Box.test(r.s^2, lag = 20, type = "Ljung-Box")
    ans = rbind(ans, c(box10[1], box10[3]))
    ans = rbind(ans, c(box15[1], box15[3]))
    ans = rbind(ans, c(box20[1], box20[3]))
    # Ljung-Box Tests - tslag required
    lag.n = 12
    x.s = as.matrix(r.s)^2
    n = nrow(x.s)
    tmp.x = .tslagGarch(x.s[, 1], 1:lag.n)
    tmp.y = x.s[(lag.n + 1):n, 1]
    fit = lm(tmp.y ~ tmp.x)
    stat = (n-lag.n) * summary.lm(fit)$r.squared
    ans = rbind(ans, c(stat, p.value = 1 - pchisq(stat, lag.n)) )
    # Add Names:
    rownames(ans) = c(
        " Jarque-Bera Test   R    Chi^2 ",
        " Shapiro-Wilk Test  R    W     ",
        " Ljung-Box Test     R    Q(10) ",
        " Ljung-Box Test     R    Q(15) ",
        " Ljung-Box Test     R    Q(20) ",
        " Ljung-Box Test     R^2  Q(10) ",
        " Ljung-Box Test     R^2  Q(15) ",
        " Ljung-Box Test     R^2  Q(20) ",
        " LM Arch Test       R    TR^2  ")
    colnames(ans) = c("Statistic", "p-Value")
    print(ans)
                               Statistic p-Value  
Jarque-Bera Test   R     Chi^2  3.021414  0.2207539
Shapiro-Wilk Test  R     W       0.9368604 0.2828191
Ljung-Box Test      R     Q(10)  8.92331     0.5393969
Ljung-Box Test      R     Q(15)  20.06886  0.1693071
Ljung-Box Test      R     Q(20)  NA             NA      
Ljung-Box Test      R^2  Q(10)  2.390212  0.9923808
Ljung-Box Test      R^2  Q(15)  6.469093  0.970746
Ljung-Box Test      R^2  Q(20)  NA             NA      
LM Arch Test         R      TR^2   5               0.957979

epoch是很好的老师，非常感谢！真的是很好的老师，就是不知道他现在在哪里高就。
本文来自: 人大经济论坛 S-Plus&R专版 版，详细出处参考：http://www.pinggu.org/bbs/viewthread.php?tid=923194&page=4&from^^uid=841068 35# zhangtao

万分感谢，真想拜你为师！ 36# epoh

好好学习一下。我也遇到上面的类似问题。郁闷。

你好，epoh。我想请问一下，当结果里出现variable：external 的时候，说明什么。这是我的Zt函数设置是t

还有一个问题，在R中，model建立好之后运行summary(lstar)错误提示为：
Error in summary.lstar(mod.lstar) :
  dims [product 224] do not match the length of object [0]
这要怎么解决啊。多谢了。

我来帮你回答下吧，不知道对不对啊，这样的提示很正常，说明你的转换变量是外部变量！希望大家有问题提出来共同解决啊！ 40# flyflyflyfly

我运行是很好啊！
41# flyflyflyfly

请问前辈，我也碰到了这样的问题
nonlinearTest(x, method="STAR-LM", p=3, d=1:2)
Problem in solve.qr(a, b): apparently singular matrix
Use traceback() to see the call stack

而且我换成其它的数据时，也出现了这样的问题，而且我把x 转换成了时间序列，问题还是存在。前辈说是数据的问题，请问有什么办法设置吗？我个人觉得不是数据的问题，我觉得是参数设置的问题，可能在nonlinearTest(x, method="STAR-LM", p=3, d=1:2)少了一个参数q吧，但是我不知道q如何设置。可能我这种想法是不对的，但是前辈能告诉我怎么解决吗？
谢谢！

数据本身的问题，中间计算用的矩阵是奇异矩阵！ 44# huaczhao

7# epoh
你好，我用了您上传的estar.r程序有个问题想请教，为什么summary之后没有t值和R方呢，能不能请您完善一下，谢谢

Epoh朋友真是高手中的高手，佩服！

真热闹，对呀，用包的时候第一次都要library一下

epoh和xuelida老师都是绝世高手，希望多指导，感谢和佩服！

请问一下你是用S-PLUS软件做的吧？在R中不能用nonlinearTest吧？ 44# huaczhao

前辈您好！我在R中估计STAR模型时出现了下面的提示：
> xx<-star(x,m=3,d=2,control=list(maxit=3000))
Using default threshold variable: thDelay=0
Testing linearity...   p-Value =  6.281377e-08
The series is nonlinear. Incremental building procedure:
Building a 2 regime STAR.
Using default threshold variable: thDelay=0
Performing grid search for starting values...
Starting values fixed: gamma =  40 , th =  10.80275 ; SSE =  576.5637
Optimization algorithm converged
Optimized values fixed for regime 2  : gamma =  40.00044 , th =  10.76126
Testing for addition of regime 3.
  Estimating gradient matrix...
  Done. Computing the test statistic...
  Done. Regime 3 is needed (p-Value =  0.01466867 ).
Adding regime  3 .
Fixing good starting values for regime  3 ...
错误于if (cost <= bestCost) { : 需要TRUE/FALSE值的地方不可以用缺少值
此外: 警告信息：
In rbind(object$model.specific$phi1, rnorm(3)) :
  number of columns of result is not a multiple of vector length (arg 2)
这是不是因为数据本身的原因出现的呢？如果不是的话怎么修正呢？谢谢！
2# epoh

Nonlinear autoregressive time series models in R
    [url]http://cran.r-project.org/web/packages/tsDyn/index.html[/url]
m: embedding dimension
d: the time delay  
s: the forecasting steps

thDelay: stands for 'delta',
         and must be an integer number between 0 and m-1.
         see page 7.(3.2 SETAR models)

mTh: stands for 'beta',
     and takes the form of a vector of real coefficients of length m.
     see page 7.(3.2 SETAR models)

mL,mH:
If not specified, mL and mH defaults to m. One can decide also only to select a
few values between 1:mL and 1:mH. This is possible using mL and mH. Hence to
have a first regime with lag 1 and 3, the second with all 3, would be: mL=c(1,3)
and mH=1:3.

phi1,phi2: see page 7.(3.2 SETAR models) formula(4)
            or page 9.(3.3 LSTAR models) generalization of SETAR

Q2:
  错误于if (cost <= bestCost)....
  方便的话上传数据,以利判断.

十分感谢，看到你在线真的很高兴，问题急需解决！下面是程序，附件里是数据。
> getwd()
[1] "d:/我的文档"
> sqr <- read.table("sqr.txt", header = TRUE);
> x <-sqr$sqr
> library(tsDyn)
>
> mod1.star <- star(x, m=4,d=1,thDelay=0,noRegimes=3,sig=0.1,control=list(maxit=3000))
Testing linearity...   p-Value =  1.665672e-05
The series is nonlinear. Incremental building procedure:
Building a 2 regime STAR.
Performing grid search for starting values...
Starting values fixed: gamma =  16 , th =  12.40319 ; SSE =  188.9531
Optimization algorithm converged
Optimized values fixed for regime 2  : gamma =  16.00069 , th =  12.39935
Testing for addition of regime 3.
  Estimating gradient matrix...
  Done. Computing the test statistic...
  Done. Regime 3 is needed (p-Value =  0.07807366 ).
Adding regime  3 .
Fixing good starting values for regime  3 ...
错误于if (cost <= bestCost) { : 需要TRUE/FALSE值的地方不可以用缺少值
此外: 警告信息：
In rbind(object$model.specific$phi1, rnorm(3)) :
  number of columns of result is not a multiple of vector length (arg 2)
>


52# epoh

还有我也遇到了下面这个问题，程序如下，数据在附件里：
> getwd()
[1] "d:/我的文档"
> library(tsDyn)
> svpdx <- read.table("quartc.txt", header = TRUE);
> x <- svpdx$qcpi
> mod1.star <- star(x, m=3,d=1,thDelay=0,noRegimes=3,sig=0.1,control=list(maxit=3000))
Testing linearity...   p-Value =  0.001151816
The series is nonlinear. Incremental building procedure:
Building a 2 regime STAR.
Performing grid search for starting values...
Starting values fixed: gamma =  21 , th =  8.19391 ; SSE =  128.2701
Optimization algorithm converged
Optimized values fixed for regime 2  : gamma =  21.00119 , th =  8.223794
Testing for addition of regime 3.
  Estimating gradient matrix...
  Done. Computing the test statistic...
  Done. Regime 3 is needed (p-Value =  0.01342169 ).
Adding regime  3 .
Fixing good starting values for regime  3 ...
Reordering regimes...
Estimating parameters of regime 3 ...
Optimized values fixed for regime  3 : gamma =  30.76407 , th =  16.5718
Optimization algorithm converged
Optimized linear values:
0.3238569 1.282448 -0.2267579 -0.1405304
-2.391707 1.009585 -1.465678 0.715712
17.98733 -1.548293 1.862583 -1.268843
Ok.
Testing for addition of regime  4 .
  Estimating gradient matrix...
  Computing the test statistic...
Regime  4  is needed (p-Value =  0.006431805 ).
Finished building a MRSTAR with  3  regimes
警告信息：
In rbind(object$model.specific$phi1, rnorm(3)) :
  number of columns of result is not a multiple of vector length (arg 2)
>

52# epoh

谢谢您啊，函数参数的问题我已经解决，我在估计下面的方程时同时出现了前面53楼54楼的错误和警告信息，数据同54楼附件的数据，程序和结果如下：
> library(tsDyn)
> svpdx <- read.table("quartc.txt", header = TRUE);
> x <- svpdx$qcpi
> mod1.star <- star(x, m=4,d=1,thDelay=0,noRegimes=3,sig=0.1,control=list(maxit=3000))
Testing linearity...   p-Value =  0.0003500272
The series is nonlinear. Incremental building procedure:
Building a 2 regime STAR.
Performing grid search for starting values...
Starting values fixed: gamma =  18 , th =  8.22092 ; SSE =  112.1446
Optimization algorithm converged
Optimized values fixed for regime 2  : gamma =  18.00000 , th =  8.24219
Testing for addition of regime 3.
  Estimating gradient matrix...
  Done. Computing the test statistic...
  Done. Regime 3 is needed (p-Value =  0.02941768 ).
Adding regime  3 .
Fixing good starting values for regime  3 ...
错误于if (cost <= bestCost) { : 需要TRUE/FALSE值的地方不可以用缺少值
此外: 警告信息：
In rbind(object$model.specific$phi1, rnorm(3)) :
  number of columns of result is not a multiple of vector length (arg 2)
>

52# epoh

修改star.R之function startingValues()
   将原来maxTh <- quantile(as.ts(s_t), .9) # percentil 90 of s_t
   改为  maxTh <- quantile(as.ts(s_t), .75)# percentil 75 of s_t
然后re-build package 'tsDyn'

########
library(tsDyn)
####data sqr
sqr <- read.table("sqr.txt", header = TRUE);
x <-sqr$sqr
mod1.star <- star(x, m=4,d=1,thDelay=0,noRegimes=3,sig=0.1,control=list(maxit=3000))

Testing linearity...   p-Value =  1.665672e-05
The series is nonlinear. Incremental building procedure:
Building a 2 regime STAR.
Performing grid search for starting values...
Starting values fixed: gamma =  16 , th =  12.40319 ; SSE =  188.9531
Optimization algorithm converged
Optimized values fixed for regime 2  : gamma =  16.00069 , th =  12.39935

Testing for addition of regime 3.
  Estimating gradient matrix...
  Done. Computing the test statistic...
  Done. Regime 3 is needed (p-Value =  0.07807366 ).
Adding regime  3 .
Fixing good starting values for regime 3 ...
Reordering regimes...
Estimating parameters of regime 3 ...
Optimized values fixed for regime  3 : gamma =  40.4925 , th =  15.80372
Optimization algorithm converged
Optimized linear values:
0.2746266 1.298432 -0.2259236 -0.04944461 -0.0386298
-2.956231 1.312587 -1.672697 0.4749626 0.1970925
23.57368 -2.076103 1.995213 -0.6978957 -0.7024029
Ok.
Testing for addition of regime  4 .
  Estimating gradient matrix...
  Computing the test statistic...
Regime  4  is NOT accepted (p-Value =  0.1812764 ).

Finished building a MRSTAR with  3  regimes
#################
####data svpdx
svpdx <- read.table("quartc.txt", header = TRUE);
xx <- svpdx$qcpi
mod2.star <- star(xx, m=4,d=1,thDelay=0,noRegimes=3,sig=0.1,control=list(maxit=3000))

Testing linearity...   p-Value =  0.0003500272
The series is nonlinear. Incremental building procedure:
Building a 2 regime STAR.
Performing grid search for starting values...
Starting values fixed: gamma =  18 , th =  8.22092 ; SSE =  112.1446
Optimization algorithm converged
Optimized values fixed for regime 2  : gamma =  18.00000 , th =  8.24219

Testing for addition of regime 3.
  Estimating gradient matrix...
  Done. Computing the test statistic...
  Done. Regime 3 is needed (p-Value =  0.02941768 ).
Adding regime  3 .
Fixing good starting values for regime 3 ...
Reordering regimes...
Estimating parameters of regime 3 ...
Optimized values fixed for regime  3 : gamma =  0.002031225 , th =  2.555149
Optimization algorithm converged
Optimized linear values:
4646656 4672.685 142.4123 -121.2756 36.74252
-66627.73 1708.731 4584.483 -841.5228 -3582.858
-9317491 120.1544 -285.2847 242.0331 -73.32033
Ok.
Testing for addition of regime  4 .
  Estimating gradient matrix...
  Computing the test statistic...
Regime  4  is NOT accepted (p-Value =  0.9999442 ).

Finished building a MRSTAR with  3  regimes

万分感谢，为什么我找不到这个文件呢？“修改star.R之function startingValues()”中的star.R，我在我的整个电脑里搜他都搜不到，谢谢您啊！
本文来自: 人大经济论坛 S-Plus&R专版 版，详细出处参考：[url=http://www.pinggu.org/bbs/viewthread.php?tid=923194&page=6&from^^uid=841068http://www.pinggu.org/bbs/viewthread.php?tid=923194&page=6&from^^uid=841068[b[/url]] 56# epoh

底下稍加说明function star()运行过程
虽然noRegimes你设置3,function star()
还是由Build the 2-regime star开始:

object <- star.predefined(x, m=m, noRegimes=2, thVar=thVar,
                         d=d, steps=steps, mTh=mTh, thDelay=thDelay,
                         series=series, trace=trace, control=control)

由于noRegimes=3,所以经过检试之后,会告诉你

Regime 3 is needed or Regime 3 is NOT accepted

如果Regime 3 is needed,就继续执行下列指令

   object <- addRegime(object);

   ("Fixing good starting values for regime")
   object <- startingValues(object, control=control, trace=trace);

   ('Estimating parameters of regime')
   object <- estimateParams(object, control=control, trace=trace);

你的问题出现在startingValues()

   startingValues.star <- function(object, trace=TRUE, ...){

   .......

   # Maximum and minimum values for c

     minTh <- quantile(as.ts(s_t), .1) # percentil 10 of s_t

     maxTh <- quantile(as.ts(s_t), .9) # percentil 90 of s_t

   .......}

maxTh太大造成cost产生NA,于是告诉你

错误于if(cost<= bestCost){:需要TRUE/FALSE值的地方不可以用缺少值

这里你可以配合你的数据调整maxTh.

star.R在 tsDyn\R
   http://cran.r-project.org/web/packages/tsDyn/index.html
         Package source:tsDyn_0.7-40.tar.gz


build package 可参考:
   http://www.murdoch-sutherland.com/Rtools/index.html


我把修改后,re-build package tsDyn_0.7-40.zip先提供你使用
   
   

tsDyn_0.7-40.zip
大小:(1.32 MB)

 马上下载

谢谢您啊！ 30# epoh

前辈您好，我又碰到新问题了，麻烦您指教。如果我选择的AR（p）模型中自回归的项数不是从1到p而是间断性的，比如说自回归中只有AR(1)和AR(3)项的话我该怎么写STAR函数呢？这样说可能不明白，举个例子，如果我想估计Yt=a+a1Y（t-1）+a2Y(t-3)+[b+b1Y（t-1）+b2Y(t-3)]G(.),其中G（.）为logistic function，如果写mod=star(x, m=3,d=1,thDelay=2,noRegimes=3,control=list(maxit=3000))的话估计出来的会把AR(2)项钱的系数估计出来，如果我只想估计AR(1)和AR(3)项的话我该如何写star（）函数呢？谢谢您了！

2# epoh