x=read.table("y.txt")
x=x[,1]
setar(x,m=2,d=2,steps=2,thDelay=0)
Non linear autoregressive model
SETAR model ( 2 regimes)
Coefficients:
Low regime:
     phiL.1      phiL.2     const L
0.08142965 -0.06675904  5.99903274
High regime:
     phiH.1      phiH.2     const H
0.03171358  0.08306024 -5.78033422
Threshold:
-Variable: Z(t) = + (1) X(t)+ (0)X(t-1)
-Value: 15.67
Proportion of points in low regime: 58.67%       High regime: 41.33%

真是太谢谢了epoh老师了，后续遇到问题再向老师请假啊，多谢指点啊

利用R data lynx,在s-plus运行ESTAR MODEL
可得结果如下:
#################in s-plus
> lynx.estar$parameters

本文来自: 人大经济论坛 S-Plus&R专版 版，详细出处参考： https://bbs.pinggu.org/forum.php? ... 1&from^^uid=11232

epoh老师，您好！
为什么我运行时提示以下错误：> lynx.estar$parameters
Problem: Object "lynx.estar" not found
Use traceback() to see the call stack
>

zhangtao兄

这一段是s-plus code,请配合看

modeling financial time series　with s-plus
CHAP 18 Nonlinear Time Series Models

#s-plus code
ESTAR.res = function(theta, g.scale, x, y, q)
{
k = ncol(x)
G = 1 - exp( - exp(theta[1])/g.scale * (q - theta[2])^2)
X = cbind(x * (1 - G), x * G)
m = crossprod(t(backsolve(chol(crossprod(X)), diag(2 * k))))
beta = m %*% t(X) %*% y
y - X %*% beta
}
data(lynx)
lynx.LHS = log10(lynx)[3:length(lynx)]
lynx.RHS = cbind(1, tslag(log10(lynx), 1:2, trim=T))
lynx.estar = nlregb(length(lynx)-2,
start=c(0,mean(lynx.RHS[,2])),
residuals=ESTAR.res,
lower=c(-Inf, min(lynx.RHS[,2])),
upper=c( Inf, max(lynx.RHS[,2])),
g.scale=var(lynx.RHS[,2]),
x=lynx.RHS, y=lynx.LHS, q=lynx.RHS[,2])
lynx.estar$parameters
exp(lynx.estar$parameters[1])/var(lynx.RHS[,2])

lynx.estar$parameters
#[1] -0.4641629 3.1881329
exp(lynx.estar$parameters[1])/var(lynx.RHS[,2])
#[1] 2.015791

+ y - X %*% beta，epoh老师，您好！以上这一句是什么意思？
为什么不是：+ y =X %*% beta？它和上一句beta = m %*% t(X) %*% y是什么关系？
nlregb这个函数如何知道它的调用方法？是splus中那个包中的函数？
另外，函数中theta没有赋值，为什么？
非常感谢！






> ESTAR.res = function(theta, g.scale, x, y, q)
+ {
+ k = ncol(x)
+ G = 1 - exp( - exp(theta[1])/g.scale * (q - theta[2])^2)
+ X = cbind(x * (1 - G), x * G)
+ m = crossprod(t(backsolve(chol(crossprod(X)), diag(2 * k))))
+ beta = m %*% t(X) %*% y
+ y - X %*% beta
+ }
> data(lynx)
> lynx.LHS = log10(lynx)[3:length(lynx)]
> lynx.RHS = cbind(1, tslag(log10(lynx), 1:2, trim=T))
> lynx.estar = nlregb(length(lynx)-2,
+ start=c(0,mean(lynx.RHS[,2])),
+ residuals=ESTAR.res,
+ lower=c(-Inf, min(lynx.RHS[,2])),
+ upper=c( Inf, max(lynx.RHS[,2])),
+ g.scale=var(lynx.RHS[,2]),
+ x=lynx.RHS, y=lynx.LHS, q=lynx.RHS[,2])
> lynx.estar$parameters
[1] -0.4641629  3.1881329
> exp(lynx.estar$parameters[1])/var(lynx.RHS[,2])
[1] 2.015791
>
> lynx.estar$parameters
[1] -0.4641629  3.1881329
> #[1] -0.4641629 3.1881329
> exp(lynx.estar$parameters[1])/var(lynx.RHS[,2])
[1] 2.015791
> #[1] 2.015791
>

nlregb这个函数如何知道它的调用方法?

就跟R一样请键入:

?nlregb
Nonlinear Least Squares Subject to Box Constraints  
    Local minimizer for sums of squares of nonlinear functions
    subject to bound-constrained parameters.

nlregb(nres, start, residuals, jacobian = NULL, scale = NULL,  
       control = NULL, lower = -Inf, upper = Inf, ...)

residuals function 就是程序中定义的ESTAR.res()

whose sum of squares is to be minimized

##########
data(lynx)
lynx.LHS = log10(lynx)[3:length(lynx)]
lynx.RHS = cbind(1, tslag(log10(lynx), 1:2, trim=T))
g.scale=var(lynx.RHS[,2])  #[1] 0.3118682
x=lynx.RHS
y=lynx.LHS
q=lynx.RHS[,2]
k = ncol(x)

##为方便说明theta就先取初始值

##theta[1] : gamma
##theta[2] : threshold
theta=c(0,mean(lynx.RHS[,2]))
#exponential transition function
G = 1 - exp( - exp(theta[1])/g.scale * (q - theta[2])^2)
X = cbind(x * (1 - G), x * G)
m = crossprod(t(backsolve(chol(crossprod(X)), diag(2 * k))))

#beta就是回归方程系数

beta = m %*% t(X) %*% y
beta
#           [,1]
#[1,]  1.5137057
#[2,]  1.4795807
#[3,] -0.9927024
#[4,]  0.7488970
#[5,]  1.1184420
#[6,] -0.3844286

#y - X %*% beta 就是residuals

y - X %*% beta

##回归方程系数beta,也可如下求出

lm.fit(X,y)

Coefficients:
              lag1       lag2              lag1       lag2
1.513706 1.479581 -0.9927024 0.748897 1.118442 -0.3844286

epoh老师，好久不见了，最近还好吧？

前辈，您好！看到您和“南冰”的帖子回复，我深深敬佩您在STAR模型和R软件方面的造诣，我现在在看这方面的东西，有很多困惑就想请教您。您做过STAR模型中脉冲响应函数的R软件实现过程吗？

你先看一下79楼

内有参考文献:

"Impulse response analysis in nonlinear multivariate models".pdf

Generalized Impulse Response Function

程序及图形.

您好!我对Ｌｓｔａｒ里面的输入参数，有些不解。。。lstar(x,m,d=1,steps=d,series,mL,mH,mTh,thDelay

m,d,steps ：embeddingdimension,timedelay,forecastingsteps　中embeddingdimension指的是什么？

前辈您好！我想请教一下您，在STAR模型中做脉冲响应分析是不是用的bootstrap模拟方法呢？我怎么有点看不懂您在讨论中给出来的那个程序，bootstrap方法不是用来求置信区间的吗？可否详细给我说说，谢谢您啦！

请先看两篇文献

看完有概念之后我好解释.

第一篇:先了解[TI] & [GI] 之区别

Smooth Transition Autoregressive Models -A Survey of Recent Developments.pdf
page 27/56
  traditional impulse response function [TI]  (38)
page 28/56
  Generalized Impulse Response Function [GI]  (41)

In case of the STAR model, analytic expressions for the conditional
expectations involved in the GI are not available for h > 1.
Stochastic simulation has to be used to obtain estimates of the
impulse response measures. See Koop et al. (1996) for a detailed
description of the relevant techniques.

接着再看第二篇的详细步骤:

Impulse response analysis in nonlinear multivariate models.pdf
Gary Koop , M. Hashem Pesaran, Simon M. Potter
5. Monte Carlo techniques for computation of GI
   page 17-19/29

epoch你好，请问R有实现LSTVAR的软件包吗？

请先参考底下两个函数

function TVAR()

   Multivariate Treshold Autoregressive model

function TVECM()

   Treshold Vector Error Correction model (VECM)

to:雁茗轩

请务必参考page 17-18/29

Monte Carlo techniques for computation of GI step1-7

同时执行79楼程序,才容易理解.

##########
N=60       #the maximum horizon
R=1000     #the number of replications
noh:112    #no. of histories
1. Pick a history and shock, W(t-1),v(t),by some combination of..
2. For a given horizon N, randomly sample (N + 1)x R values
   of the (Kdimensional)innovation
   #generate standardized shocks from lstar model
    st_shock_idx=sample(seq(1:length(std_res)),R*(N+1),replace = TRUE)   
    length(st_shock_idx)
    st_shock=matrix(std_res[st_shock_idx],R,N+1)
    dim(st_shock)  # 1000 x 61
    shockz=st_shock*sd_res
3. Use the first N random shocks (obtained under step 2) to
   compute the realization.....
4. Use the same draw of N random shocks plus one additional
   draw of the random shock to produce a realization...
5. Repeat steps 3 and 4 R times and form the averages...
   ###
   for(i in 1:noh){
   hist=hist_m[i,];
   # benchmark profile   
   .....
   .....
   ######compute "new" histories
   histv=cbind(1,y,histv[,2],y,histv[,2])
   realvb[k]=mean(histv[,4]);
   }  #end k shock profile
6. Take the difference between the two averages to form a
   Monte Carlo estimate of the GI,
   GI[i,]=realvb-realzb;
   }# end i numbers of histories

太感谢您了，我会尽快看完再请教您，谢谢！

epoh老师您好，我想请问一下您，在R软件中估计出STAR模型的各个参数之后，怎么能得到线性部分的系数的标准差呢？谢谢您了！

epoh老师您好，我想请问一下您，在R软件中估计出STAR模型的各个参数之后，怎么能得到线性部分的系数的标准差呢？谢谢您了！

73楼

我想请问下thDelay 到底是什么意思呢？

非常好的帖子。学习。

能不能作双变量的LSTAR回归呢。

同样的数据，下面二个命令结果有什么不同？？
1、mod.lstar<-lstar(x,m=2,d=2,control=list(maxit=3000))
Testing linearity...   p-Value =  0.2626276
The series is linear. Using linear model instead.

2、mod<-star(x,m=2,d=2,thDelay=1,noRegimes=3,control=list(maxit=3000))
Using maximum autoregressive order for low regime: mL = 2
Using maximum autoregressive order for high regime: mH = 2
Using default threshold variable: thDelay=0
Performing grid search for starting values...
Starting values fixed: gamma =  40 , th =  2.5762 ; SSE =  562.6991
Optimization algorithm converged
Optimized values fixed for regime 2  : gamma =  57.80452 , th =  2.563136

为什么出现不同的结果，能为我说明么？谢谢！

不对呢：
你说到：“其实你也可以使用function star
因为star 也是使用Logistic transition function
底下两者结果相同
mod.star <- star(log10(lynx), m=2, d=1, control=list(maxit=3000))
mod.star
summary(mod.star)

mod.lstar <- lstar(log10(lynx), m=2, d=1, control=list(maxit=3000))
mod.lstar”
      但我用同样的数据(Y)进行对比，出现不同的结果；如下。
     > mod.star <- star(y, m=2, d=1, control=list(maxit=3000))
Using default threshold variable: thDelay=0
Testing linearity...   p-Value =  0.4377235
The series is linear. Using linear model instead.
      > mod.star
Non linear autoregressive model
AR model
Coefficients:
       const        phi.1        phi.2
0.001580686  1.213115336 -0.327246860
   
     > mod.lstar <- lstar(y, m=2, d=1, control=list(maxit=3000))
Using maximum autoregressive order for low regime: mL = 2
Using maximum autoregressive order for high regime: mH = 2
Using default threshold variable: thDelay=0
Performing grid search for starting values...
Starting values fixed: gamma =  40 , th =  0.9588955 ; SSE =  43.58875
Optimization algorithm converged
Optimized values fixed for regime 2  : gamma =  41.8667 , th =  0.9560731
    > mod.lstar
   Non linear autoregressive model
   LSTAR model
Coefficients:
Low regime:
     phi1.0      phi1.1      phi1.2
-0.05123484  1.12120401 -0.29997985

High regime:
     phi2.0      phi2.1      phi2.2
0.45453226  0.05158612 -0.14314690
Smoothing parameter: gamma = 41.87
Threshold
Variable: Z(t) = + (1) X(t) + (0) X(t-1)
Value: 0.956
>
这可是质的差别呢。




本文来自: 人大经济论坛 S-Plus&R专版 版，详细出处参考： https://bbs.pinggu.org/forum.php? ... 3&from^^uid=57021

数据是log10(lynx)时,底下两者相同
  mod.star <- star(log10(lynx), m=2, d=1, control=list(maxit=3000))
    Using default threshold variable: thDelay=0
    Testing linearity...   p-Value =  0.001858152
    The series is nonlinear.
  mod.lstar <- lstar(log10(lynx), m=2, d=1, control=list(maxit=3000))

换上你的数据
mod.star <- star(y, m=2, d=1, control=list(maxit=3000))
Using default threshold variable: thDelay=0
Testing linearity...   p-Value =  0.4377235
The series is linear. Using linear model instead.

    在62楼中分别对thDelay取不同的值，如
    thDelay=0,  thDelay=1 这有什么讲研么？？

thDelay : 'time delay' for the threshold variable

log10(lynx):
  2.429752
  2.506505
  2.767156
  ....
  ....
  3.201397
  3.424392
  3.530968

##########
thDelay=0
thVar
  2.506505
  2.767156
  2.940018
  ....
  ....
  3.000000
  3.201397
  3.424392

##########
thDelay=1
thVar
  2.429752
  2.506505
  2.767156
  ....
  ....
  2.820858
  3.000000
  3.201397

##########
thDelay=0
  2.506505
  2.767156
  2.940018
  ....
  ....
  3.000000
  3.201397
  3.424392

你好；我还是搞不清，用我的数据会得到不同的结果！
1、用mod.star得到线性的结果。
   > mod.star <- star(y, m=2, d=1, control=list(maxit=3000))
Using default threshold variable: thDelay=0
Testing linearity...   p-Value =  0.4377235
The series is linear. Using linear model instead.
>   mod.star
Non linear autoregressive model（非线性自回归）
AR model
Coefficients:
       const        phi.1        phi.2
0.001580686  1.213115336 -0.327246860 ；
2、用mod.lstar（可以得到二机制）
    > mod.lstar <- lstar(y, m=2, d=1, control=list(maxit=3000))
Using maximum autoregressive order for low regime: mL = 2
Using maximum autoregressive order for high regime: mH = 2
Using default threshold variable: thDelay=0（用系统默认的方法来确定转换变量或转换门限值！）
Performing grid search for starting values...
Starting values fixed: gamma =  40 , th =  0.9588955 ; SSE =  43.58875
Optimization algorithm converged
Optimized values fixed for regime 2  : gamma =  41.8667 , th =  0.9560731
    > mod.lstar
Non linear autoregressive model
LSTAR model
Coefficients:
Low regime:
     phi1.0      phi1.1      phi1.2
-0.05123484  1.12120401 -0.29997985
High regime:
     phi2.0      phi2.1      phi2.2
0.45453226  0.05158612 -0.14314690
Smoothing parameter: gamma = 41.87
Threshold
Variable: Z(t) = + (1) X(t) + (0) X(t-1)
Value: 0.956
     这二者有什么差别呢，在实际我该用那一种办法？是不是说不能作STAR，而是作LSTAR呢？
3、如果M取高些，如M=5，那么在回归中会现会出回归系数不显著的情况（线性和非线性部分中回归项），如何去分别这些不显著的系数呢：如果能分别出其中的不显著系数，如phi1.3和phi2.4系数不显著，那么又何去设定方程呢！如设定为：
   mod.lstar <- lstar(y, m=6, d=1, control=list(maxit=3000))，如何发现有不显著的系数，又如何变换这个式子呢。
    4、thDelay是不是指的门限变量的滞后阶数呢，也就是说如果在多变量的LSTAR中，以其它解释变量的滞后项作为开关变量时，可通过thDelay来设定滞后阶数！如X（-5）为开关变量，则thDelay=6！
   谢谢！

    你好：“epoh”
    我真的很学会这种方法。你能不能就此开一个专题，专门讲LSTAR或STAR用R。在讲的过程中一天一讲，讲解到每个参数的设置，留出一些问题给各位复习和操作。一星期一大练习，练到活学活用；在讲的时候也用些中文，我们可理解的快一些。

我也有这样的问题，不知道如何解决。
    1、m=3，是不是控制自回归的最高阶数，即AR（P），P<=m;
    2、thDelay是不是控制转换（开关变量）变量的最大滞后阶数，如S(t-d)？
    3、回归结果的说明中，说了最优化的结果是这样的。那么，在各区中是不是就不存在ar(p)回归系数不显著的情况！通过设置了系数，系统默认这个结果（也就是没有象在EVIEWS中那样对各回归系数的显著性进行检验），也就没有对各系数显者性进行检验的必要！
    4、如果不是以上的说法（结果），那么正如你所说的：如果检验AR(1)和AR(3)项系数显著的话，那该如何写star（）函数！