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

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

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=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: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
#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')

前辈您好！您可以把您在下面这个帖子的29楼中提到的
您修改后的lstar_new.R的code发给我吗？因为我又遇到了28楼的情况，我想通过修改minTh和maxTh来实现我的估计，如果方便的话可以把.R格式的lstar（）函数和star()函数的code也发给我一份可以吗？谢谢您了！这样的话我可以通过source（）来调用函数，也可以修改minTh和maxTh，万分感谢！
https://bbs.pinggu.org/thread-923194-3-1.html
29# epoh

内含四个文件:
修改过 : lstar_new.R , isLinear.R
原始文件 : lstar.R , star.R


lstar files.rar (21.59 KB)

2011-5-21 09:02:51 上传

前辈您好！我修改maxTh时是不是也得对应的修改minTh呢？比如说以前
初始值为：  # Maximum and minimum values for c
    minTh <- quantile(as.ts(z), .1) # percentil 10 de z     
    maxTh <- quantile(as.ts(z), .9) # percentil 90 de z     
    rateTh <- (maxTh - minTh) / 200;      
如果我想将maxTh修改为0.75，是不是必须同时将minTh修改为0.25呢？即：
# Maximum and minimum values for c
    minTh <- quantile(as.ts(z), .25) # percentil 25 de z     
    maxTh <- quantile(as.ts(z), .75) # percentil 75 de z     
    rateTh <- (maxTh - minTh) / 200;      
谢谢您了！
本文来自: 人大经济论坛 S-Plus&R专版 版，详细出处参考：https://bbs.pinggu.org/viewthread.php?tid=923194&page=6&from^^uid=841068
56# epoh

你说的问题作者解决了
请下载最新版 May 15, 2011
tsDyn_0.7-60

刚看了source code
你28楼所说的问题,
作者新写了一个function calculateLinearCoefficients()来解决,
原来 #newPhi1 <- lm(yy ~ . - 1, data.frame(tmp))$coefficients;
改后 newPhi1 <- calculateLinearCoefficients(tmp, yy)

summary(mod.lstar) 也已更新,可以使用.
但isLinear()可能忘记更新了,还是用我给你的.

万分感谢啊，估计star model的问题如您所说的已经解决，但是最近我做estar model时又碰到了如下提示“错误于if (cost <= bestCost) { : 需要TRUE/FALSE值的地方不可以用缺少值”，您可以参照修改过的star（）这个函数的程序把您以前发给我的estar.R这个函数的程序修改下吗？谢谢您了啊！ 85# epoh

estar.R
estar_new.rar (5.22 KB) 本附件包括：

• estar.R

2011-5-21 20:10:38 上传

前辈您好！麻烦您检查下发给我的这个estar_new.Ri这个函数的程序可以吗？因为我用两组数据估计estar model，但是这两组数据都存在估计出的有些系数的结果非常大的现象，不大合情理，太谢谢您了！
下面是估计出的结果：
> mod=estar(y, m=2, d=1,thDelay=1,trace=TRUE, control=list(3000))
Using maximum autoregressive order for low regime: mL = 2
Using maximum autoregressive order for high regime: mH = 2
Performing grid search for starting values...
Starting values fixed: gamma =  28 , th =  1.9871 ; SSE =  80.93593
Optimization algorithm converged
Optimized values fixed for regime 2  : gamma =  28 , th =  1.9871
> mod
Non linear autoregressive model
ESTAR model
Coefficients:
Low regime:
     phi1.0      phi1.1      phi1.2
2916974.72   -44116.55 -1810883.56
High regime:
     phi2.0      phi2.1      phi2.2
-2916974.41    44117.39  1810883.14
Smoothing parameter: gamma = 28
Threshold
Variable: Z(t) = + (0) X(t) + (1) X(t-1)
Value: 1.987
>


87# epoh

重新检查一次
1.Transition function没错.
  同s-plus 公式(18.18)
   G <- function(y, g, th)  1-exp(- g*(y-th)^2)  

2.function gradEhat()需要配合修改
  这是用在line 163
  res <- optim(p, SS, gradEhat, hessian = TRUE, method="BFGS", control = control)
  但是也可设为NULL,这时optim就会自动使用finite-difference approximation
  所以 line 163 请你改为
  res <- optim(p, SS, NULL, hessian = TRUE, method="BFGS", control = control)
3.借用你的数据验证
  estar有时系数的确很大
  但是fitted values还OK
4.Results:
svpdx <- read.table("data.txt", header = TRUE);
x <- svpdx$cpi    #19
y<-svpdx$rpi
z<-svpdx$m2
##################data cpi
ndx.lstar <- lstar(x, m=3,d=1, thVar=z,control=list(maxit=3000));
ndx.estar <- estar(x, m=3,d=1, thVar=z,control=list(maxit=3000));
#####################
ndx.lstar$str$yy
t(ndx.lstar$fitted.values)
t(ndx.estar$fitted.values)

> ndx.lstar$str$yy
[1] 24.1 17.1  8.3  2.8 -0.8 -1.4  0.4  0.7 -0.8  1.2  3.9  1.8  1.5  4.8
[15]  5.9 -0.7
> t(ndx.lstar$fitted.values)
         [,1]     [,2]     [,3]     [,4]      [,5]      [,6]      [,7]
[1,] 20.90693 17.21399 8.527687 1.814728 0.3473067 0.3291664 0.7412002
         [,8]       [,9]     [,10]    [,11]    [,12]     [,13]    [,14]
[1,] 1.174782 -0.9689601 -2.808649 3.083912 4.392991 -1.426489 2.629319
        [,15]    [,16]
[1,] 7.584261 5.257824

> t(ndx.estar$fitted.values)
         [,1]     [,2]    [,3]     [,4]      [,5]      [,6]      [,7]
[1,] 20.44798 17.09968 8.30002 2.798391 0.8791174 -2.078809 0.3993102
         [,8]        [,9]    [,10]    [,11]    [,12]     [,13]    [,14]
[1,] 2.203886 -0.06964867 -2.58948 3.554523 5.220085 -1.813132 1.828305
        [,15]    [,16]
[1,] 7.685169 4.934606

##################data rpi
ndx1.lstar <- lstar(y, m=3,d=1, thVar=z,control=list(maxit=3000));
ndx1.estar <- estar(y, m=3,d=1, thVar=z,control=list(maxit=3000));
ndx1.lstar$str$yy
t(ndx1.lstar$fitted.values)
t(ndx1.estar$fitted.values)

> ndx1.lstar$str$yy
[1] 21.7 14.8  6.1  0.8 -2.6 -3.0 -1.5 -0.8 -1.3 -0.1  2.8  0.8  1.0  3.8
[15] 5.9 -1.2
> t(ndx1.lstar$fitted.values)
         [,1]     [,2]     [,3]       [,4]      [,5]      [,6]       [,7]
[1,] 18.10048 14.91605 6.472074 -0.2602880 -1.495639 -1.973182 -0.7282473
           [,8]      [,9]     [,10]    [,11]   [,12]     [,13]    [,14]
[1,] -0.7556872 -1.684395 -2.502178 1.543447 4.00788 -2.471291 2.296770
        [,15]   [,16]
[1,] 5.726296 6.00791

> t(ndx1.estar$fitted.values)
         [,1]     [,2]     [,3]      [,4]      [,5]       [,6]      [,7]
[1,] 21.68860 14.80551 5.560748 0.9267582 -2.617242 -0.6988458 -1.393348
           [,8]       [,9]      [,10]     [,11]    [,12]      [,13]    [,14]
[1,] -0.1743817 -0.1635246 -0.4696633 0.5452012 2.060080 -0.1212813 1.199337
        [,15]    [,16]
[1,] 2.714618 3.337430

前辈您好！非常感谢您我一直以来对我的帮助，还有几个问题希望得到您的解答！
1、在修改程序时，是不是只需将line 163 的“res <- optim(p, SS, gradEhat, hessian = TRUE, method="BFGS", control = control)”
   修改为“res <- optim(p, SS, NULL, hessian = TRUE, method="BFGS", control = control)"就可以吗？还是同时需要修改前面的”function gradEhat()“呢？：
gradEhat <- function(p)
    {
      gamma <- p[1]  #Extract parms from vector p
      th    <- p[2]       #Extract parms from vector p
      new_phi<- lm.fit(cbind(xxL, xxH * G(z, gamma, th)), yy)$coefficients
      phi1 <- new_phi[1:(mL+1)]
      phi2 <- new_phi[(mL+2):(mL + mH + 2)]
      y.hat <- F(phi1, phi2, gamma, th)
      e.hat <- yy - y.hat
      fX <- sigmoid(gamma * (z - th));
      dfX <- dsigmoid(fX);
      
      gGamma <- as.vector(xxH %*% phi2) * as.vector(dfX * (z - th));
      gTh <-        - as.vector(xxH %*% phi2) * as.vector(gamma * dfX);
      J = - cbind(gGamma, gTh) / sqrt(str$n.used)
      
      return(2 * t(e.hat) %*% J)
      
    }
因为我只将line 163 的“res <- optim(p, SS, gradEhat, hessian = TRUE, method="BFGS", control = control)”
   修改为“res <- optim(p, SS, NULL, hessian = TRUE, method="BFGS", control = control)"但是却得不到您估计的结果，只得到了如下结果：
> source("estar.R")
> svpdx <- read.table("data.txt", header = TRUE);
> x <- svpdx$cpi    #19
> y<-svpdx$rpi
> z<-svpdx$m2
> ndx.estar <- estar(x, m=3,d=1, thVar=z,control=list(maxit=3000));
Using maximum autoregressive order for low regime: mL = 3
Using maximum autoregressive order for high regime: mH = 3
Using only first 16 elements of thVar
Performing grid search for starting values...
Starting values fixed: gamma =  40 , th =  0.2139298 ; SSE =  118.9796
Optimization algorithm converged
Optimized values fixed for regime 2  : gamma =  299.3408 , th =  0.5413974
> ndx.estar
Non linear autoregressive model
ESTAR model
Coefficients:
Low regime:
    phi1.0     phi1.1     phi1.2     phi1.3
  554456.1 -5265009.3  4618125.2 -1483153.7
High regime:
    phi2.0     phi2.1     phi2.2     phi2.3
-554456.3  5265011.2 -4618126.6  1483154.1
Smoothing parameter: gamma = 299.3
Threshold
Variable: external
Value: 0.5414
>
2、第二个问题是在调用data.txt数据后：
svpdx <- read.table("data.txt", header = TRUE);
x <- svpdx$cpi    #19
y<-svpdx$rpi
z<-svpdx$m2
模型估计时ndx.estar <- estar(x, m=3,d=1, thVar=z,control=list(maxit=3000));转换变量z只用了前16个数据，其实正确的应该是用后面16个数据，因为x（t）=phi1.0 +phi1.1*x（t-1）+phi1.2*x（t-2）+phi1.3*x（t-3）+（phi2.0 +phi2.1*x（t-1）+phi2.2*x（t-2）+phi2.3*x（t-3））*G(z（t）)
，但是我只用后面是六个数据时系统会报错：
”错误于lm.fit(xx, yy) : 外接函数调用时不能有NA/NaN/Inf(arg1)“
3、我的新数据如下，麻烦您帮我调试下程序重新估计下，谢谢您了！
程序如下：
source("estar.R")
svpdx1 =read.table("m.txt", header = TRUE);
y=svpdx1$m2
y
svpdx2 =read.table("mzh.txt", header = TRUE);
x=svpdx2$dm
x
mod=estar(y, m=2, d=1,thVar=x, trace=TRUE, control=list(3000))；
89# epoh