我把exponential function修改为
s-plus 公式(18.22) ,加入g.scale
感觉比较稳定.
跑出的结果,几乎跟s-plus相同

gamma = 1.013 , s-plus 1.013556
th=-2.775 , s-plus -2.774638

##########################

source("estar_gscale.R")
ndx <- read.table("ndx.rvol.txt");
ndx=ndx[,1]
mod.estar <- estar(log(ndx), m=2,  control=list(maxit=3000))
mod.estar
ESTAR model
Coefficients:
Low regime:
     phi1.0      phi1.1      phi1.2
-0.02512237  0.82408380  0.19486576
High regime:
     phi2.0      phi2.1      phi2.2
-2.14572870 -0.74500270  0.05390954
Smoothing parameter: gamma = 1.013
Threshold
Variable: Z(t) = + (1) X(t) + (0) X(t-1)
Value: -2.775

estar_gscale

estar_gscale.rar
大小:(3.36 KB)

 马上下载

本附件包括：

• estar_gscale.R

Q2:
external threshold variable
程序会重第一个算起
若要改变可以自己指定thVar=z[4:19]
ndx.estar <- estar(x, m=3,d=1, thVar=z[4:19],control=list(maxit=3000));
ndx.estar$model.specific$thVar

Q3:
新数据的问题来自你的 x
x当external threshold variable是否合适??
> y
[1]  4.678046 25.864895 26.321201 21.249118 15.900649 18.735037
> x
[1]  0.4563057 -5.0720827 -5.3484693  2.8343885 16.1039869 -9.4525352

补充:
希望你上传的数据,
不是你要研究的数据,
因为都无法拒绝Linearity Tests


source("isLinear.R")
########data lynx
mod.lstar <- lstar(log10(lynx), m=2, mTh=c(0,1), control=list(maxit=3000))
isLinear(mod.lstar)   #reject
#Using default threshold variable: thDelay=0
#firstOrderTest thirdOrderTest fifthOrderTest
#   0.000365692    0.001858152    0.020038384
##########s-plus data
ndx <- read.table("ndx.rvol.txt");
ndx=ndx[,1]
mod.ndx <- lstar(log10(ndx), m=2,  control=list(maxit=3000))
isLinear(mod.ndx)     #reject
#Using default threshold variable: thDelay=0
#firstOrderTest thirdOrderTest fifthOrderTest
#   0.037127869    0.001443528    0.001939517
########data cpi
svpdx <- read.table("data.txt", header = TRUE);
x <- svpdx$cpi    #19
y<-svpdx$rpi
z<-svpdx$m2
ndx.lstar <- lstar(x, m=2,d=1, thVar=z[4:19],control=list(maxit=3000));
isLinear(ndx.lstar)
#Using default threshold variable: thDelay=0
#firstOrderTest thirdOrderTest fifthOrderTest
#    0.38286906     0.40174922     0.08667536
######## data m2 & dm
svpdx1 =read.table("m.txt", header = TRUE);
y=svpdx1$m2
y
svpdx2 =read.table("mzh.txt", header = TRUE);
x=svpdx2$dm
x
mod2=lstar(y,m=2,d=1,thVar=x,trace=TRUE,control=list(3000))
isLinear(mod2)
#Using default threshold variable: thDelay=0
#firstOrderTest thirdOrderTest fifthOrderTest
#     0.5099226      0.9538843      0.8437448

好的，万分感谢啊，不过这个模型还是可以做LSTAR模型的，isLinear检验的仅仅是
“Using default threshold variable: thDelay=0”
> library(tsDyn)
> svpdx1 =read.table("m.txt", header = TRUE);
> y=svpdx1$m2
> y
[1]  4.678046 25.864895 26.321201 21.249118 15.900649 18.735037 34.839024
[8] 25.386489 29.219643 22.849062 22.379741 18.316832 27.974895 26.528477
[15] 31.276486 37.312023 34.529817 29.466797 25.257197 19.580787 14.840376
[22] 14.736983 12.269493 13.579229 19.856236 19.634701 15.500372 17.570000
[29] 16.950000 16.700000 17.820000 27.700000
> svpdx2 =read.table("mzh.txt", header = TRUE);
> x=svpdx2$dm
> x
[1]  0.4563057 -5.0720827 -5.3484693  2.8343885 16.1039869 -9.4525352
[7]  3.8331535 -6.3705804 -0.4693215 -4.0629090  9.6580637 -1.4464183
[13]  4.7480087  6.0355369 -2.7822062 -5.0630191 -4.2096000 -5.6764102
[19] -4.7404113 -0.1033927 -2.4674896  1.3097354  6.2770067 -0.2215349
[25] -4.1343289  2.0696282 -0.6200000 -0.2500000  1.1200000  9.8800000
> mod=star(y, m=2, d=1,thVar=x, trace=TRUE, control=list(3000))
Testing linearity...   p-Value =  7.035431e-15
The series is nonlinear. Incremental building procedure:
Building a 2 regime STAR.
Performing grid search for starting values...
Starting values fixed: gamma =  40 , th =  2.836254 ; SSE =  171.8763
Optimization algorithm converged
Optimized values fixed for regime 2  : gamma =  40.0004 , th =  2.842359
Testing for addition of regime 3.
  Estimating gradient matrix...
  Done. Computing the test statistic...
  Done. Regime 3 is needed (p-Value =  3.669726e-08 ).
Adding regime  3 .
Fixing good starting values for regime  3 ...
Reordering regimes...
Estimating parameters of regime 3 ...
Optimized values fixed for regime  3 : gamma =  41.76118 , th =  6.513988
Optimization algorithm converged
Optimized linear values:
3.220247 0.832652 -0.1018466
0.5934429 0.4479709 -0.1300669
-96.12018 6.00583 -0.3571657
Ok.
Testing for addition of regime  4 .
  Estimating gradient matrix...
  Computing the test statistic...
Regime  4  is needed (p-Value =  1.668348e-25 ).
Adding regime  4 .
Fixing good starting values for regime  4 ...
Reordering regimes...
Estimating parameters of regime 4 ...
Optimized values fixed for regime  4 : gamma =  45.94539 , th =  10.64555
Optimization algorithm converged
Optimized linear values:
-1223.773 0.9989231 0.001131645
2450.004 0.002154858 -0.002266246
-0.000291771 1.323034e-05 5.945269e-06
8.071782e-07 1.512423e-05 1.283766e-05
Ok.
Testing for addition of regime  5 .
  Estimating gradient matrix...
  Computing the test statistic...
Regime  5  is needed (p-Value =  7.299121e-05 ).
Adding regime  5 .
Fixing good starting values for regime  5 ...
Reordering regimes...
Estimating parameters of regime 5 ...
Optimized values fixed for regime  5 : gamma =  40.19385 , th =  12.99012
Optimization algorithm converged
Optimized linear values:
-87770.04 0.9999943 -1.059852e-06
175536.9 1.140156e-05 2.120564e-06
8.209452e-06 4.957738e-06 -6.154125e-06
-1.467837e-06 -5.366139e-06 6.189698e-06
-3.195256e-07 -6.826985e-06 8.100021e-06
Ok.
Testing for addition of regime  6 .
  Estimating gradient matrix...
  Computing the test statistic...
Regime  6  is NOT accepted (p-Value =  0.9999734 ).
Finished building a MRSTAR with  5  regimes
>


92# epoh

前辈您好！我还想做出两个负的冲击下的的impulse response function图(这个是不是只需将“shock_m=2; #Double Magnitude Positive Shock”改为“shock_m=-2; #Double Magnitude negative Shock”就可以了呢？)和两个正的冲击与两个负的冲击之和的impulse response function图（想衡量非对称的脉冲响应Measuring asymmetric impulse response，见帖子中26楼的paper“Smooth Transition Autoregressive Models -A Survey of Recent Developments”page28），麻烦您帮我下编程序，谢谢您了！
81# epoh

estar transition functions
依 s-plus 公式(18.22)加入g.scale
感觉比较稳定
跑出的结果,几乎跟s-plus相同.
程序我放91楼.

######################
asymmetric impulse response
麻烦你告知要用的数据及mode,以方便编程
因为tsDyn_0.7-60还是有些改变,
这次shock_m我会改为,一次完成
shock_m=c(-2,-1,0,1,2);
当然你也可以改为跟文献一样
shock_m=seq(-3,3,0.1)

#####################
m,mzh的数据OK,是我忘了加thVar=x
svpdx1 =read.table("m.txt", header = TRUE);
y=svpdx1$m2
svpdx2 =read.table("mzh.txt", header = TRUE);
x=svpdx2$dm
mod2=lstar(y,m=2,d=1,thVar=x,trace=TRUE,control=list(3000))
isLinear(mod2,thVar=x)    #reject
firstOrderTest thirdOrderTest fifthOrderTest
2.041730e-211  3.925313e-155  3.347176e-101

非常感谢！数据和模型都和您在81楼的相同，只是用现在的tsDyn_0.7-60估计出来的结果变了，我还是想用以前的这个模型和结果(如下所示)做asymmetric impulse response function。谢谢您了！
library(tsDyn)
svpdx <- read.table("quartc.txt", header = TRUE);
x <- svpdx$qcpi
mod<- star(x, m=5,d=1,thDelay=1,noRegimes=3, control=list(maxit=3000))
Testing linearity...   p-Value =  0.002852162
The series is nonlinear. Incremental building procedure:
Building a 2 regime STAR.
Performing grid search for starting values...
Starting values fixed: gamma =  24 , th =  14.00072 ; SSE =  88.80465
Optimization algorithm converged
Optimized values fixed for regime 2  : gamma =  23.99999 , th =  13.99887
Testing for addition of regime 3.
  Estimating gradient matrix...
  Done. Computing the test statistic...
  Done. Regime 3 is needed (p-Value =  0.01449991 ).
Adding regime  3 .
Fixing good starting values for regime 3 ...
Reordering regimes...
Estimating parameters of regime 3 ...
Optimized values fixed for regime  3 : gamma =  42.27714 , th =  13.93399
Optimization algorithm converged
Optimized linear values:
0.2997076 1.056419 0.06167327 0.04290917 -0.5146651 0.2624537
1.415723 1.049638 -1.271572 -0.5466812 0.785233 -0.1815548
12.89779 -1.462572 1.594120 0.3835062 -1.037747 -0.005988798
Ok.
Testing for addition of regime  4 .
  Estimating gradient matrix...
  Computing the test statistic...
Regime  4  is NOT accepted (p-Value =  0.1873368 ).
Finished building a MRSTAR with  3  regimes
> mod
Non linear autoregressive model
Multiple regime STAR model
Regime  1 :
    Linear parameters: 0.2997076, 1.0564187, 0.0616733, 0.0429092, -0.5146651, 0.2624537
Regime  2 :
    Linear parameters: 1.4157232, 1.0496377, -1.2715717, -0.5466812, 0.785233, -0.1815548
    Non-linear parameters:
35.0696739, 5.2725434
Regime  3 :
    Linear parameters: 12.8977894, -1.4625725, 1.5941205, 0.3835062, -1.0377474, -0.0059888
    Non-linear parameters:
42.2771382, 13.9339929
95# epoh

library(tsDyn)
svpdx =read.table("quartc.txt", header = TRUE);
x =svpdx$qcpi
mod =star(x, m=5,d=1,thDelay=1,noRegimes=3,control=list(maxit=3000))
phi1=mod$model.specific$phi1
phi2=mod$model.specific$phi2
res=mod$residuals;
sd_res=sd(res)
std_res=res/sd_res;
hist_m=cbind(1,mod$str$xx,mod$str$yy,mod$model.specific$thVar)
colnames(hist_m)=c("Const","V1/0","V1/-1","V1/-2","V1/-3",",V1/-4","y","thVar")
hist_m[1:3,]         #"histories"
noh=dim(hist_m)[1]   # 83 "histories"
#Transition function
   G <- function(y, g, th)  plogis(y, th, 1/g)
N=60       #the maximum horizon
R=1000     #the number of replications
shock_m=c(-2,-1,0,1,2);
nosh=length(shock_m)
realvb=NA*seq(1:(N+1));
realzb=NA*seq(1:(N+1));
GI=matrix(data = NA,noh,N+1);
GIRF= array(0, dim = c(noh,N+1, nosh))
#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)
st_shock[1:2,]
for(s in 1:nosh){
shockz=st_shock*sd_res
###
for(i in 1:noh){
hist=hist_m[i,];
# benchmark profile
histz=matrix(rep(hist,R),R,dim(hist_m)[2], byrow =TRUE);
dim(histz)
histz[1:3,]
for(j in 1:(N+1)){
y=histz[,1:6]%*%phi1[1,] + histz[,1:6]%*%phi1[2,]*G(histz[,8],phi2[1,1],phi2[1,2]) +
  histz[,1:6]%*%phi1[3,]*G(histz[,8], phi2[2,1],phi2[2,2])+shockz[,j]
######compute "new" histories
histz=cbind(1,y,histz[,2:5],y,histz[,3])
realzb[j]=mean(histz[,7]);
}#end j benchmark profile

#shock profile
shockv=matrix(rep(shock_m*sd_res,R),R,1)
histv=matrix(rep(hist,R),R,dim(hist_m)[2], byrow =TRUE)
k=1
y=histz[,1:6]%*%phi1[1,] + histz[,1:6]%*%phi1[2,]*G(histz[,8],phi2[1,1],phi2[1,2]) +
  histz[,1:6]%*%phi1[3,]*G(histz[,8], phi2[2,1],phi2[2,2])+shockv
######compute "new" histories
histv=cbind(1,y,histv[,2:5],y,histv[,3])
realvb[k]=mean(histv[,7]);
for(k in 2:(N+1)){
........
........
}  #end k shock profile
GI[i,]=realvb-realzb;
}# end i numbers of histories
GIRF[,,s]=GI
} #end s
girf1=apply( GIRF[,,1],2,mean)
girf2=apply( GIRF[,,2],2,mean)
girf3=apply( GIRF[,,3],2,mean)
girf4=apply( GIRF[,,4],2,mean)
girf5=apply( GIRF[,,5],2,mean)
par(mfrow = c(2, 2))

##P shock
plot(girf4, ylab = expression(GIRF(h,delta,W(t-1))), xlab = 'Time Horizon')
title("Generalized impulse response functions (P shock)")
legend("bottomright", " Positive Shock",text.col=4)
lines(girf4,col='red')

##PP shock
plot(girf5, ylab = expression(GIRF(h,delta,W(t-1))), xlab = 'Time Horizon')
title("Generalized impulse response functions (PP shock)")
legend("bottomright", "Double Magnitude Positive Shock",text.col=4)
lines(girf5,col='red')

##N shock
plot(girf2, ylab = expression(GIRF(h,delta,W(t-1))), xlab = 'Time Horizon')
title("Generalized impulse response functions (N shock)")
legend("bottomright", " Negative Shock",text.col=4)
lines(girf2,col='red')

##NN shock
plot(girf1, ylab = expression(GIRF(h,delta,W(t-1))), xlab = 'Time Horizon')
title("Generalized impulse response functions (NN shock)")
legend("bottomright", "Double Magnitude Negative Shock",text.col=4)
lines(girf1,col='red')




####Measure of asymmetry.
x11()
asym=girf1+girf5  #NN + PP
plot(asym, ylab = expression(ASYM(h,delta,W(t-1))), xlab = 'Time Horizon')
title("Measure of asymmetry")
legend("bottomright", "Double Magnitude Positive Negative Shock",text.col=4)
lines(asym,col='red')