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

我还是以data lynx来做图形
正负冲击对比鲜明,
这样比较贴近你的想法.
###########PP shock
gi_pphigh=GIRF[,,5][which(hist_m[,4]>th),]    #34
gi_pplow =GIRF[,,5][which(hist_m[,4]<th),]    #78
##high pp shock
high_regime5=apply(gi_pphigh ,2,mean)
##low pp shock
low_regime5=apply(gi_pplow ,2,mean)
####plot PP shock on different regime
par(mfrow = c(2,1))
plot(seq(1,61),high_regime5, type="l", col = 2, ylab = expression(GIRF(h,delta,W(t-1))), xlab = 'Time Horizon')
lines(seq(1,61),low_regime5,col = 3, lty = 2)
title("PP shock on different regimes")
legend("bottomright", c("High Regime","Low Regime"), lty=1:2, col=2:3,    adj = c(0, .6))
###############
###########NN shock
gi_nnhigh=GIRF[,,1][which(hist_m[,4]>th),]    #34
gi_nnlow =GIRF[,,1][which(hist_m[,4]<th),]    #78
##low NN shock
##high nn shock
high_regime1=apply(gi_nnhigh ,2,mean)
##low nn shock
low_regime1=apply(gi_nnlow ,2,mean)
####plot NN shock on different regimes
plot(seq(1,61),high_regime1, type="l", col = 2,ylab = expression(GIRF(h,delta,W(t-1))), xlab = 'Time Horizon')
lines(seq(1,61),low_regime1,col = 3, lty = 2)
title("NN shock on different regimes")
legend("bottomright", c("High Regime","Low Regime"), lty=1:2, col=2:3,    adj = c(0, .6))



###############plot PP &NN Shock on high regime
x11()
par(mfrow = c(2,1))
plot(seq(1,61),high_regime5, type="l", col = 2,ylim=c(-0.6,0.6),
   ylab = expression(GIRF(h,delta,W(t-1))), xlab = 'Time Horizon')
lines(seq(1,61),high_regime1,col = 3, lty = 2)
title("PP & NN shock on high regime")
legend("bottomright", c("PP Shock","NN Shock"), lty=1:2, col=2:3,    adj = c(0, .6))
###############plot PP &NN Shock on low regime
plot(seq(1,61),low_regime5, type="l", col = 2,  ylim=c(-0.55,0.55),
    ylab = expression(GIRF(h,delta,W(t-1))), xlab = 'Time Horizon')
lines(seq(1,61),low_regime1,col = 3, lty = 2)
title("PP & NN shock on low regime")
legend("bottomright", c("PP Shock","NN Shock"), lty=1:2, col=2:3,    adj = c(0, .6))




ps:你的数据三区的话
gi_pphigh=GIRF[,,5][which(hist_m[,7]>phi2[2,2]),]                    #9
gi_ppmiddle=GIRF[,,5][which(hist_m[,7]<phi2[2,2] & hist_m[,7]>phi2[1,2]),]  #15
gi_pplow=GIRF[,,5][which(hist_m[,7]<phi2[1,2]),]   

前辈您好！可以把你的例子数据data lynx传给我吗？谢谢您了！101# epoh

抱歉没说清楚
这是延续79楼
data(lynx)就是package范例
我把程序贴上好了.

##################
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=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: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)){
........
........
}  #end k shock profile
GI[i,]=realvb-realzb;

}# end i numbers of histories
GIRF[,,s]=GI
} #end s

续接101楼

前辈您好！我现在安装的是R 2.13.0，直接估计下面方程会得到如下结果：
> 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.99876
Testing for addition of regime 3.
  Estimating gradient matrix...
  Done. Computing the test statistic...
  Done. Regime 3 is NOT accepted (p-Value =  0.06172199 ).
Finished building a MRSTAR with  2  regimes
> phi1=mod$model.specific$phi1
> phi1
           [,1]       [,2]       [,3]        [,4]       [,5]        [,6]
[1,]  0.1883403  1.4341040 -0.3197945 -0.13549328 -0.1913063  0.16429109
[2,] 14.4249331 -0.7906222  0.7040176  0.01522778 -0.5758736 -0.08938245
但是我还是想用R 2.12.2估计模型以得到下面的结果，
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
但是我在R 2.12.2中选择在线安装的话还是得到第一种结果，如果选择本地安装tsDyn_0.7-40.zip，但是估计时出现了下面的问题：
> utils:::menuInstallLocal()
程序包'tsDyn'打开成功，MD5和检查也通过
> library(tsDyn)
错误: 程辑包'tsDyn'没有安装在'arch=i386'中
>  svpdx <- read.table("quartc.txt", header = TRUE);
>  x <- svpdx$qcpi
>  mod2 <- star(x, m=5,d=1,thDelay=1,noRegimes=3,control=list(maxit=3000))
错误: 没有"star"这个函数
>
现在我还是想估计得到第二种结果，该怎么办呢？谢谢您了！
103# epoh

安装 2.13.0 也是OK
我就是安装2.13.0
你把现有的R移除乾净(含package)
然后安装2.13.0
由本机安装03.23我上传的版本
载入tsDyn,他会告诉你缺甚么,
你就安装什么.
如此就OK了.


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

你想画highest density regions for generalized impulse responses
可以安装Package "hdrcde"
Author 就是Rob J Hyndman
底下就是我用data(lynx)画的

epoh真是好老师，向您学习！

前辈您好！可以把您做的highest density regions for generalized impulse responses的程序贴上来吗？谢谢您了啊！
105# epoh

hdr.boxplot_modify.R

hdr.boxplot_modify.rar (1.75 KB) 本附件包括：

• hdr.boxplot_modify.R

2011-6-2 08:38:10 上传

epoh您好，我在使用tsDyn程序包里面的TVAR.LRtest做两个还是三个门限检测的时候出现了问题说错误于object$CriticalValBoot[1, ] : 量度数目不对，能否麻烦您帮我看一下怎么回事啊，谢谢了啊！
> data(zeroyld)
> data<-zeroyld
> test<-TVAR.LRtest(data,lag=2,mTh=1,thDelay=1:2,nboot=3,plot=FALSE,trim=0.1,test="2vs3")
Warning: the thDelay values do not correspond to the univariate implementation in tsdyn
> summary(test)
Test of linear AR against TAR(1) and TAR(2)

LR test:
          [,1]
Test  12.69123
P-Val  1.00000
Bootstrap critical values for test 1 vs 2 regimes
错误于object$CriticalValBoot[1, ] : 量度数目不对

epoh您好，tsDyn里面有没有对TVECM是2个门限和3个门限检验的检验程序啊？如果我想使用Ming Chien Lo和Eric Zivot(2001)的方法进行2VS3的门限检验，请问怎么可以得到检验量近视分布的统计量和P值啊？多谢多谢！