在“Nonlinear autoregressive time series models in R usingtsDynversion 0.7”  ast revised 11/03/2008 by Antonio, Fabio Di Narz中，我用系统自带的例子和数据。The first model proposed in the literature for these data was an AR(2)。
采用：
mod.ar <- linear(x, m=2)
mod.ar
summary(mod.ar)
Coefficient(s):
       Estimate  Std. Error  t value  Pr(>|t|)   
const  1.057600    0.120837   8.7523 2.957e-14 ***
phi.1  1.384238    0.063332  21.8569 < 2.2e-16 ***
phi.2 -0.747776    0.063385 -11.7973 < 2.2e-16 ***
可以看出，T值检验都过；但是如果我不是M=2，而是直接设M=9等。也可通过T值来检验显著性。如；
Coefficient(s):
       Estimate  Std. Error  t value  Pr(>|t|)   
const  0.761169    0.288317   2.6400  0.009661 **
phi.1  1.224610    0.095749  12.7898 < 2.2e-16 ***
phi.2 -0.665402    0.152797  -4.3548 3.305e-05 ***
phi.3  0.275708    0.165320   1.6677  0.098595 .  
phi.4 -0.328840    0.166685  -1.9728  0.051363 .  
phi.5  0.158732    0.166583   0.9529  0.343023   
phi.6 -0.130592    0.165276  -0.7901  0.431370   
phi.7  0.071094    0.153231   0.4640  0.643711   
phi.8  0.131305    0.096537   1.3601  0.176937   
也就是说有些不显著，那么我修改参数M=4
> mod.ar <- linear(x, m=4)
> summary(mod.ar)
Coefficient(s):
       Estimate  Std. Error  t value  Pr(>|t|)   
const  1.429164    0.190764   7.4918 2.229e-11 ***
phi.1  1.270993    0.093915  13.5335 < 2.2e-16 ***
phi.2 -0.702846    0.153078  -4.5914 1.226e-05 ***
phi.3  0.146591    0.153094   0.9575   0.34050   
phi.4 -0.206564    0.094052  -2.1963   0.03027 *  
还有PHI。3不显著，那么要如何改参数M呢？我试了一下，设M=C（1、2，4）可不行呢。如下，
> mod.ar <- linear(x, m=c(1,2,4))
错误于structure(c(this, list(...)), class = class(this)) : 下标出界
此外: 警告信息：
1: In if ((m * d + steps) > n) stop("time series too small to handle these embedding parameters") :
  条件的长度大于一，因此只能用其第一元素
2: In 0:(m - 1) : 数值表达式一共有3元素: 只用了第一个
3: In 1:m : 数值表达式一共有3元素: 只用了第一个
> 下标出界是什么意思！该如何修改M参数！谢谢！

现样：在上文中；我运用autopairsn也报错（第13页）。
〉autopairs(x, lag=1, type="regression", GUI=FALSE)
Loading required package: sm
错误于autopairs(x, lag = 1, type = "regression", GUI = FALSE) :
  sm package is required for kernel estimations
此外: 警告信息：
In library(package, lib.loc = lib.loc, character.only = TRUE, logical.return = TRUE,  :
  there is no package called 'sm'

Q : 我试了一下，设M=C（1、2，4）可不行呢?
      mod.ar <- linear(x, m=c(1,2,4))
ANS: m : embedding dimension
               is a scalar NOT a vector

Q : 如果检验AR(1)和AR(3)项系数显着的话，那该如何写star（）函数！
ANS: 需要依自己模型,自行编程

Q : autopairs(x, lag=1, type="regression", GUI=FALSE)
ANS: function autopairs() need package "sm"
        pls install "sm" package

对一个在职的人是很难的

顶一个  回头细看

请问各位前辈，在R软件中估计STAR模型时，引入外生变量做转换变量怎么估计呢？比如，现在x是被解释变量，引入y做转换变量，怎么写代码呢？很着急啊，请各位慷慨赐教，谢谢啦！mod=star(x,m=2,d=1,thDelay=1,noRegimes=2,control=list(maxit=3000))

哈哈!参考一下28楼(data)及36楼

老师您好，我就是按照那样估计的，可是为什么mod=lstar(x,m=2,d=1,thvar=z,thDelay=0和mod=lstar(x,m=2,d=1,thDelay=0）估计出来的模型一样，查看了系数和thVar值，都是x的数据，那z压根没有起作用吗？

我想请教你个问题，关于STAR模型的估计的 ，y是被解释变量，我想引入外生变量x作为转换变量进行估计，可是得到的结果与用y的滞后期作为转换变量的结果、以及x的差分值作为转换变量的结果一样，这是怎么回事呢？
比如mod=star(y,m=3,d=1,thDelay=0,thVar=x,noRegimes=2)与mod=star(y,m=3,d=1,thDelay=0,noRegimes=2)结果是一样的？麻烦问一下这是怎么回事呢？谢谢前辈了！

这是数据和具体的程序，麻烦您看一下问题出在哪里，CPI是被解释变量Y,Oil是解释变量x。

老师，数据和程序在160楼，麻烦您看一下我的问题在哪里，万分感谢了，着急用啊！

library(tsDyn)
data=read.table("data.txt", header = TRUE)
x <- data$cpi
y <- data$rpi
z <- data$m2
mod=lstar(x,m=2,d=1,control=list(maxit=3000))
mod

Non linear autoregressive model

LSTAR model
Coefficients:
Low regime:
    const1     phi1.1     phi1.2
1.3866542  0.5369498 -0.4856777

High regime:
    const2     phi2.1     phi2.2
15.6621561  0.2037306 -0.5327252

Smoothing parameter: gamma = 22.35

Threshold
Variable: Z(t) = + (1) X(t) + (0) X(t-1)

Value: 6.348

mod1=lstar(x,m=3,d=1,thVar=z,control=list(maxit=3000))
mod1

Non linear autoregressive model

LSTAR model
Coefficients:
Low regime:
   const1    phi1.1    phi1.2    phi1.3
-0.739900  1.852660 -1.256096  0.731050

High regime:
    const2     phi2.1     phi2.2     phi2.3
2.7868699 -1.0567380  1.0397138 -0.8612963

Smoothing parameter: gamma = 330

Threshold
Variable: external
Value: 0.2821

mod1$model.specific$thVar
[1]  0.26524514  0.31278198  0.37310154  0.34529154  0.29467111 -0.03867458
[7]  0.34668927  0.17245702  0.14411363 -0.31433668  0.87572994  0.16337324
[13]  0.19679602  0.18161695  0.14387323  0.18886459

library(tsDyn)
nan=read.table("cpo.txt", header = TRUE)
y=nan$CP
x=nan$Oil
mod=lstar(y,m=3,d=1,thVar=x,control=list(maxit=3000))
mod

Non linear autoregressive model

LSTAR model
Coefficients:
Low regime:
      const1       phi1.1       phi1.2       phi1.3
0.213101468  0.311218447 -0.147596767  0.003313981

High regime:
    const2     phi2.1     phi2.2     phi2.3
-0.2364208 -0.3765884  1.7958783  0.3482814

Smoothing parameter: gamma = 34

Threshold
Variable: external
Value: 93.04

mod$model.specific$thVar
  [1]  19.72  20.72  24.53  26.18  27.04  25.52  26.97  28.39  29.66  28.84

真是太感谢您了，epoh老师！！可是为什么我的就是估计不出来呢?难道是R软件中安装的有问题？还有，老师可否把程序或者估计出来的残差贴给我呢，谢谢啦！还有，我看帖子前面您和别的坛友交流时可以计算参数t值和标准差等，为什么我的也出来这些东西呢？谢谢啦！

最重要的是,要在你的电脑能运行,才能进行研究工作
我的版本是
R version 2.15.0
tsDyn Version: 0.8-1

你先执行看看有问题再说

好的，我重新安装一下试试……

请教epoh,
Splus 与R包tsDyn中的LSTAR 模型估计出来的 gamma差异在哪里？
其它参数差不多，但gamma似乎算的不是一个gamma。
比如lynx数据
Splus：      STAR(log10(lynx),p=2,d=2)
R：         lstar(log10(lynx),m=2,thDelay=1)      

貌似 tsDyn直接给出的是 1/(1+exp(-gamma * (z-th))) 中的gamma
Splus 给出的是 1/(1+exp(-gamma/sd(z) *(z-th))) 中的gamma

呵呵,iooo兄是行家.
的确是这样.
Transition function

  R： 1/(1+exp(-gamma * (z-th)))
  Splus : 1/(1+exp(-gamma/sd(z) *(z-th)))
  Splus 多了个transition.scale  gamsc=stdev(z)

不过我个人觉得R比较正确
R是用grid search 来寻找gamma
这点可由fitted.values看得出来
个人观点仅供参考

library(tsDyn)
nan=read.table("cpo.txt", header = TRUE)
y=nan$CP
x=nan$Oil
mod=lstar(y,m=2,d=1,thVar=x,control=list(maxit=3000))
mod
xxL=cbind(1,mod$str$xx)
xxH=xxL
yy=mod$str$yy
z=mod$model.specific$thVar
gamma=mod$coefficients[7]
th=mod$coefficients[8]

#Transition function
    G <- function(y, g, th)  plogis(y, th, 1/g)

#lm.fit(cbind(xxL, xxH * G(z, gamma, th)), yy)

calc.lm.t <- function(x)
{
   Qr <- x$qr
   r <- x$residuals
   p <- x$rank
   p1 <- 1L:p
   rss <- sum(r^2)

   n <- NROW(Qr$qr)
   rdf <- n - p

   resvar <- rss/rdf
   R <- chol2inv(Qr$qr[p1, p1, drop = FALSE])
   se <- sqrt(diag(R) * resvar)

   est <- x$coefficients[Qr$pivot[p1]]
   tval <- est/se

   res <- cbind(est = est, se = se, tval = tval)
   res
}
lmf <- lm.fit( cbind(xxL, xxH * G(z, gamma, th)), yy)
calc.lm.t(lmf)

#                  est            se                 tval
#           0.02619398 0.09426718  0.2778696
#V1/0   0.30976606 0.13298950  2.3292521
#V1/-1 -0.13333347 0.13577256 -0.9820355
#          0.25543047 0.12957109  1.9713539
#V1/0  -0.07042287 0.18376493 -0.3832226
#V1/-1  0.10601972 0.18798123  0.5639910

mod
Non linear autoregressive model

LSTAR model
Coefficients:
Low regime:
     const1      phi1.1      phi1.2
0.02619145  0.30976991 -0.13333290

High regime:
     const2      phi2.1      phi2.2
0.25543391 -0.07042918  0.10601783

Smoothing parameter: gamma = 12.42

Threshold
Variable: external
Value: 56.38

恩，谢谢。没有怎么仔细的去看代码

真的非常感谢您，您是德艺双馨的前辈！我重新安装了一遍R和程序包，可以估计了

老师您好！我把模型估计出来了，想对残差做个ARCH-LM检验，检验时，总是出现要建立一个varest之类的东西，这个要怎么解决呢，谢谢了！老麻烦您真不好意思……

library(tsDyn)
nan=read.table("cpo.txt", header = TRUE)
y=nan$CP
x=nan$Oil
mod=lstar(y,m=2,d=1,thVar=x,control=list(maxit=3000))
mod
names(mod)
#[1] "str"      "coefficients"   "fitted.values"  "residuals"   "k"            
#[6] "model"    "model.specific"

#########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 = mod$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  2.378766  0.3044091
Shapiro-Wilk Test  R    W      0.9891328 0.4601097
Ljung-Box Test     R    Q(10)  6.890205  0.735768  
Ljung-Box Test     R    Q(15)  25.69184  0.0413892
Ljung-Box Test     R    Q(20)  33.93066  0.02659757
Ljung-Box Test     R^2  Q(10)  9.998199  0.4406513
Ljung-Box Test     R^2  Q(15)  10.39486  0.7941945
Ljung-Box Test     R^2  Q(20)  13.69091  0.8458153
LM Arch Test       R    TR^2   13.27055  0.3496864

非常感谢您，我好好看看！您真是个好老师……

老师您好，为什么我自己做的时候，参照您给出来的程序写好，总是出现下面的错误呢：
错误于`rownames<-`(`*tmp*`, value = c(" Jarque-Bera Test   R    Chi^2 ",  :
  'dimnames'的长度[1]必需与陈列范围相等
>      colnames(ans) = c("Statistic", "p-Value")
>      print(ans)
     Statistic p-Value   
[1,] 0.6940057 0.7068033  
[2,] 0.9950685 0.9540148  
[3,] 7.852276  0.6432638  
[4,] 32.20009  0.00604615
[5,] 38.24426  0.008266291
[6,] 5.550793  0.8514875  
[7,] 10.21307  0.8061273  
[8,] 14.72797  0.7917557
而且列名并不显示是哪个统计量？应该出来9个统计量的，这只有8个，是不是我哪个地方开始设置错了？
麻烦您了！

程序你有修改吧?
请把修改的程序贴出来或贴在短信息

数据是我自己估计出来的一个模型，和您上面估计的那个只是换了一个变量而已，样本容量一致，然后在下面的估计中我把残差命名为res,您命名的是r.s.其他的都没有改变，然后出来的结果就是这样子……

> #########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")

Residuals Tests:
>      res = mod$residuals
>      ans = NULL
>  # Normality Tests:
>      jbtest = jarqueberaTest(res)@test
>      ans = rbind(ans, c(jbtest[1], jbtest[2]))
>      if (length(res) < 5000) {
+          swtest = shapiro.test(res)
+          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(res, lag = 10, type = "Ljung-Box")
>      box15 = Box.test(res, lag = 15, type = "Ljung-Box")
>      box20 = Box.test(res, 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(res^2, lag = 10, type = "Ljung-Box")
>      box15 = Box.test(res^2, lag = 15, type = "Ljung-Box")
>      box20 = Box.test(res^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(res)^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)$resquared
>      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  ")
错误于`rownames<-`(`*tmp*`, value = c(" Jarque-Bera Test   R    Chi^2 ",  :
  'dimnames'的长度[1]必需与陈列范围相等
>      colnames(ans) = c("Statistic", "p-Value")
>      print(ans)
     Statistic p-Value   
[1,] 0.6940057 0.7068033  
[2,] 0.9950685 0.9540148  
[3,] 7.852276  0.6432638  
[4,] 32.20009  0.00604615
[5,] 38.24426  0.008266291
[6,] 5.550793  0.8514875  
[7,] 10.21307  0.8061273  
[8,] 14.72797  0.7917557

nan_res.txt
   
   

真是不好意思，太麻烦你了老师！