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

前辈您好！再请教你个问题啊。我boot三区制LSTAR模型的code如下：
#by bootstrapping residuals
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[1,]
phi1
phi2=mod$model.specific$phi1[2,]
phi2
phi3=mod$model.specific$phi1[3,]
phi3
re=mod$residuals
mre=(re-mean(re))/sd(re)
n=length(x)    #88
n1=length(mre)  #83
number=0*seq(1:n1)
B=3
gamma1=c(1:B)
gamma2=c(1:B)
#Transition function
   G1 <- function(y1, g1, th1)  plogis(y1, th1, 1/(g1))
   G2 <- function(y2, g2, th2)  plogis(y2, th2, 1/(g2))

for  (b in 1:B) {
#bootstrap by resampling the residuals
#sampling with replacement from the original sample
  rand=sample(seq(1:n1),n1, replace = TRUE)
  rr=mre[rand]
#generate new y
  y=c(x[1:5], 6 :n)
      for  (t in 6:n){
  
   y[t]=t(c(1,y[(t-1):(t-5)]))%*%phi1[1,] + (t(c(1,y[(t-1):(t-5)]))%*%phi1[2,])*G(y[t-2],
phi2[1,1],phi2[1,2]) + (t(c(1,y[(t-1):(t-5)]))%*%phi1[3,])*G(y[t-2], phi2[2,1],phi2[2,2])+rr[t
-5]}
  
#{y[t]=0.2997076+1.0564187*y[t-1]+0.0616733*y[t-2]+0.0429092*y[t-3]-0.5146651*y[t-4]+0.2624537*y[t-5]+(1.4157232+1.0496377*y[t-1]-1.2715717*y[t-2]-0.5466812*y[t-3]+0.785233*y[t-4]-0.1815548*y[t-5])/(1+exp(-35.0696739*(y[t-2]-5.2725434)))+(12.8977894-1.4625725*y[t-1]+1.5941205*y[t-2]+0.3835062*y[t-3]-1.0377474*y[t-4]-0.0059888*y[t-5])/(1+exp(-42.2771382*( y[t-2]-13.9339929)))+rr[t]}#generate new
y
z=y
rs=star(z, m=5,d=1,thDelay=1,noRegimes=3,trace=FALSE,control=list(maxit=3000))
if (is.null(rs$coefficients[19]) == TRUE) gamma1=NA  else
      gamma1=rs$coefficients[19]
if (is.null(rs$coefficients[20]) == TRUE) gamma2=NA  else
      gamma2=rs$coefficients[20]
}
gamma1=gamma1[!is.na(gamma1)]          # remove NAs
gamma2=gamma2[!is.na(gamma2)]          # remove NAs
length(gamma1)     
ength(gamma2)     
但是给了我错误提示：
> #by bootstrapping residuals
>
> 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
>
> phi1=mod$model.specific$phi1[1,]
> phi1
[1]  0.29970762  1.05641871  0.06167327  0.04290917 -0.51466513  0.26245373
>
> phi2=mod$model.specific$phi1[2,]
> phi2
[1]  1.4157232  1.0496377 -1.2715717 -0.5466812  0.7852330 -0.1815548
>
> phi3=mod$model.specific$phi1[3,]
> phi3
[1] 12.897789438 -1.462572488  1.594120467  0.383506197 -1.037747375
[6] -0.005988798
>
> re=mod$residuals
> mre=(re-mean(re))/sd(re)
> n=length(x)    #88
> n1=length(mre)  #83
> number=0*seq(1:n1)
> B=3
> gamma1=c(1:B)
> gamma2=c(1:B)
>
> #Transition function
>    G1 <- function(y1, g1, th1)  plogis(y1, th1, 1/(g1))
>    G2 <- function(y2, g2, th2)  plogis(y2, th2, 1/(g2))
>
>
> for  (b in 1:B) {
+ #bootstrap by resampling the residuals
+ #sampling with replacement from the original sample
+   rand=sample(seq(1:n1),n1, replace = TRUE)
+   rr=mre[rand]
+ #generate new y
+   y=c(x[1:5], 6 :n)
+       for  (t in 6:n){
+   
+    y[t]=t(c(1,y[(t-1):(t-5)]))%*%phi1[1,] + (t(c(1,y[(t-1):(t-5)]))%*%phi1[2,])*G(y[t-2],
+ phi2[1,1],phi2[1,2]) + (t(c(1,y[(t-1):(t-5)]))%*%phi1[3,])*G(y[t-2], phi2[2,1],phi2[2,2])+rr[t
+ -5]}
+   
+ #{y[t]=0.2997076+1.0564187*y[t-1]+0.0616733*y[t-2]+0.0429092*y[t-3]-0.5146651*y[t-4]+0.2624537*y[t-5]+(1.4157232+1.0496377*y[t-1]-1.2715717*y[t-2]-0.5466812*y[t-3]+0.785233*y[t-4]-0.1815548*y[t-5])/(1+exp(-35.0696739*(y[t-2]-5.2725434)))+(12.8977894-1.4625725*y[t-1]+1.5941205*y[t-2]+0.3835062*y[t-3]-1.0377474*y[t-4]-0.0059888*y[t-5])/(1+exp(-42.2771382*( y[t-2]-13.9339929)))+rr[t]}#generate new
+ y
+ z=y
+ rs=star(z, m=5,d=1,thDelay=1,noRegimes=3,trace=FALSE,sig=1,control=list(maxit=3000))
+ if (is.null(rs$coefficients[19]) == TRUE) gamma1=NA  else
+       gamma1=rs$coefficients[19]
+
+ if (is.null(rs$coefficients[20]) == TRUE) gamma2=NA  else
+       gamma2=rs$coefficients[20]
+ }
错误于phi1[1, ] : 量度数目不对
> gamma1=gamma1[!is.na(gamma1)]          # remove NAs
> gamma2=gamma2[!is.na(gamma2)]          # remove NAs
> length(gamma1)     # 66
[1] 3
> ength(gamma2)     # 66
麻烦您帮我检查下code，谢谢您了！顺便和您说下，你上次发给我的程序有点小问题，不过我改了后可以运行了，万分感谢啊！
70# 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
phi1
#> phi1
#           [,1]      [,2]        [,3]        [,4]       [,5]         [,6]
#[1,]  0.2997076  1.056419  0.06167327  0.04290917 -0.5146651  0.262453729
#[2,]  1.4157232  1.049638 -1.27157169 -0.54668116  0.7852330 -0.181554807
#[3,] 12.8977894 -1.462572  1.59412047  0.38350620 -1.0377474 -0.005988798

phi2=mod$model.specific$phi2
phi2
#        gamma        th
#[1,] 35.06967  5.272543
#[2,] 42.27714 13.933993


re=mod$residuals
mre=(re-mean(re))/sd(re)
n=length(x)    #88
n1=length(mre)  #83
number=0*seq(1:n1)
B=10
gamma1=c(1:B)
gamma2=c(1:B)
#Transition function
   G <- function(y, g, th)  plogis(y, th, 1/(g))

for  (b in 1:B) {
#bootstrap by resampling the residuals
#sampling with replacement from the original sample
  rand=sample(seq(1:n1),n1, replace = TRUE)
  rr=mre[rand]
#generate new y
  y=c(x[1:5], 6 :n)
      for  (t in 6:n){
  
   y[t]=t(c(1,y[(t-1):(t-5)]))%*%phi1[1,] + (t(c(1,y[(t-1):(t-5)]))%*%phi1[2,])*G(y[t-2],
phi2[1,1],phi2[1,2]) + (t(c(1,y[(t-1):(t-5)]))%*%phi1[3,])*G(y[t-2], phi2[2,1],phi2[2,2])+rr[t
-5]}
  
z=y
rs=star(z, m=5,d=1,thDelay=1,noRegimes=3,control=list(maxit=3000))
if (is.null(rs$coefficients[19]) == TRUE) gamma1=NA  else
      gamma1=rs$coefficients[19]
if (is.null(rs$coefficients[20]) == TRUE) gamma2=NA  else
      gamma2=rs$coefficients[20]
}
gamma1=gamma1[!is.na(gamma1)]          # remove NAs
gamma2=gamma2[!is.na(gamma2)]          # remove NAs
length(gamma1)     
length(gamma2)     

程序修改后,可以运行,不过不稳.
th很容超出范围,造成NA.

借个角落,来答覆,耽搁许久,
短信息询问系数标准差问题的朋友nmmd

library(tsDyn)
mod =star(log10(lynx), m=2,d=1,noRegimes=2,control=list(maxit=3000))
xxL=cbind(1,mod$str$xx)
xxH=xxL
yy=mod$str$yy
z=mod$model.specific$thVar
gamma=mod$model.specific$phi2[1]
th=mod$model.specific$phi2[2]

#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.4055717 0.3096907  1.309603
#V1/0   1.2378257 0.1548691  7.992723
#V1/-1 -0.3265076 0.1040744 -3.137252
#       0.8187712 0.3604585  2.271472
#V1/0   0.3087342 0.1743599  1.770672
#V1/-1 -0.6418812 0.1319384 -4.865009
mod

Non linear autoregressive model
Multiple regime STAR model
Regime  1 :
    Linear parameters: 0.4055717, 1.2378257, -0.3265076

Regime  2 :
    Linear parameters: 0.8187712, 0.3087342, -0.6418812

    Non-linear parameters:
40.0002197, 2.5661053

前辈，您好！我现在是想做STAR的Impulse response function，但找不到impulse（）这个函数
，但是在google中搜到了附件中的这个文件（见我的帖子的跟帖），作者做STAR model时是调用
的RSTAR这个函数包，但是tsDyn这个包中没有impulse（）、NlStarTest（）等函数，而这些函数
在附件中提到的RSTAR包中却有。Smooth Transition Autoregressive Models - A
Survey Of Recent development这篇文章在这个帖子的26楼，这篇文章的page26也有Impulse response function的
介绍，麻烦您帮我看下怎么在R中实现实现LSTAR模型的Impulse response function，谢谢您了啊!  
73# epoh

##########
page 13/56
Testing linearity against LSTAR
有function isLinear()可用
mod.lstar <- lstar(log10(lynx), m=2, mTh=c(0,1), control=list(maxit=3000))
mod.lstar
isLinear(mod.lstar)

#Using default threshold variable: thDelay=0
#firstOrderTest thirdOrderTest fifthOrderTest
#   0.000365692    0.001858152    0.020038384

75# epoh

epho老师，我是听tulipsliu老师说起您，来围观了下您的帖子，真的好佩服啊！学术和为人都是很多人不可及！向您学习！

能在论坛遇到高人指教也很幸运，羡慕各位鞋童！

祝福每一位热心人！

AZA AZA Fighting!

学习了 多谢楼主啊

epoh老师绝对是大牛啊，人品学术没的说啊，谢谢epoh老师啊！ 76# O(∩_∩)O~！

Generalized Impulse Response Function
  page 18/29
  "Impulse response analysis in nonlinear
  multivariate models".pdf
data:lynx (114 obs)
     山猫（Lynx）与雪兔（Snowshoe Hare）
     掠食与猎物,9~10年周期性消长
N=60      #the time horizon
R=1000    #the number of replications
shock_m=2;#Double Magnitude Positive Shock
histories:112

###########
###########
library(tsDyn)
mod <- lstar(log10(lynx), m=2, mTh=c(0,1), control=list(maxit=3000))
phi1=mod$model.specific$coefficients[1:3]
phi1
phi2=mod$model.specific$coefficients[4:6]
phi2
gamma=mod$model.specific$coefficients[7]
th=mod$model.specific$coefficients[8]
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","y","thVar")
hist_m[1:3,]   #"histories"
noh=dim(hist_m)[1]
#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=2; #Double Magnitude Positive Shock
nosh=length(shock_m)
realvb=NA*seq(1:(N+1));
realzb=NA*seq(1:(N+1));
GI=matrix(data = NA,noh,N+1);
#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,]
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:3]%*%phi1 + histz[,1:3]%*%phi2*G(histz[,5], gamma,th) + shockz[,j]
######compute "new" histories
histz=cbind(1,y,histz[,2],y,histz[,2])
realzb[j]=mean(histz[,4]);
}#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=histv[,1:3]%*%phi1 + histv[,1:3]%*%phi2*G(histv[,5], gamma,th) + shockv
######compute "new" histories
histv=cbind(1,y,histv[,2],y,histv[,2])
realvb[k]=mean(histv[,4]);
for(k in 2:(N+1)){
y=histv[,1:3]%*%phi1 + histv[,1:3]%*%phi2*G(histv[,5], gamma,th) + shockz[,k]
######compute "new" histories
histv=cbind(1,y,histv[,2],y,histv[,2])
realvb[k]=mean(histv[,4]);
}  #end k shock profile
GI[i,]=realvb-realzb;
}# end i numbers of histories
#GI
girf=apply(GI,2,mean)
plot(girf, ylab = expression(GIRF(h,delta,W(t-1))), xlab = 'Time Horizon')
title("Generalized impulse response functions (PP shock)")
legend("topright", "Double Magnitude Positive Shock",text.col=4)
lines(girf,col='red')

前辈您好，非常感谢您给我编的程序，不过还得麻烦您，您写的程序有些地方我看不大懂，所以我不知道三区制AR(5)的impulse function的code如何在发您给我的code的基础上修改，STAR model的code如下：
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))


数据还是54楼的数据，麻烦您帮我改下程序可以吗？谢谢您了！

79# epoh