急：R软件加载程序包tsDyn问题

我把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

2011-5-24 11:37:51 上传

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')

前辈您好！我观察了脉冲响应函数图，有一点不是很明白。就是给系统一个正的冲击的脉冲响应函数图和给系统一个负的冲击的脉冲响应函数图在零时刻时应该大小相同正负相反，以后的时刻可以不对称（也正是体现了不对称性），但是图中的情况并不是这样，希望前辈解答，谢谢！paper“Smooth Transition Autoregressive Models -A Survey of Recent Developments”page55）中的图我贴在了98楼已作对照！
97# epoh

这两种图形是不一样的
附图是highest density regions using density quantile method

要分清楚当然是依
recession and expansion分别画出

由于你是分三区:
   high  :9
  middle : 15
   low   : 59

  high regime几乎不受正负冲击影响
  为何会这样?这要依据你的数据了

我画出middle regime 给你参考

前辈您好！十分感谢您一直以来对我的帮助，您可以再把高区制及低区制的impulse response function也帖出来吗？另外把三个区制的impulse response function的code也发给我好吗？我想再看些文献分析下，谢谢您了！
99# epoh