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

Note:this modified package only for
         function star() m=4,ar(1),ar(3)
tsDyn_0.7-40.zip (1.32 MB)

2011-3-29 19:46:21 上传

library(tsDyn)
####data sqr
sqr <- read.table("sqr.txt", header = TRUE);
x <-sqr$sqr
mod1.star <- star(x, m=4,d=1,thDelay=0,noRegimes=3,sig=0.1,control=list(maxit=3000))
[1] "original"
     V1/0 V1/-1 V1/-2 V1/-3
[1,] 2.14  1.33  1.31  1.14
[2,] 1.70  2.14  1.33  1.31
[1] "ar(1),ar(3)"
     V1/0 V1/-1 V1/-3
[1,] 2.14  1.33  1.14
[2,] 1.70  2.14  1.31
Testing linearity...   p-Value =  1.665672e-05
The series is nonlinear. Incremental building procedure:
Building a 2 regime STAR.
Performing grid search for starting values...
Starting values fixed: gamma =  10 , th =  12.48815 ; SSE =  192.9971
Optimization algorithm converged
Optimized values fixed for regime 2  : gamma =  9.978157 , th =  12.46917
Testing for addition of regime 3.
  Estimating gradient matrix...
  Done. Computing the test statistic...
  Done. Regime 3 is NOT accepted (p-Value =  0.4722626 ).
Finished building a MRSTAR with  2  regimes
> mod1.star
Non linear autoregressive model
Multiple regime STAR model
Regime  1 :
    Linear parameters: 0.2333271, 1.5366889, -0.5899859, 0.0212266
Regime  2 :
    Linear parameters: 12.4843681, -0.5445538, 0.2590598, -0.4629664
    Non-linear parameters:
9.9781565, 12.4691698

#####
####data svpdx
svpdx <- read.table("quartc.txt", header = TRUE);
xx <- svpdx$qcpi
mod2.star <- star(xx, m=4,d=1,thDelay=0,noRegimes=3,sig=0.1,control=list(maxit=3000))
[1] "original"
     V1/0 V1/-1 V1/-2 V1/-3
[1,] 7.53 15.40 24.73 27.63
[2,] 4.03  7.53 15.40 24.73
[1] "ar(1),ar(3)"
     V1/0 V1/-1 V1/-3
[1,] 7.53 15.40 27.63
[2,] 4.03  7.53 24.73
Testing linearity...   p-Value =  0.0003500272
The series is nonlinear. Incremental building procedure:
Building a 2 regime STAR.
Performing grid search for starting values...
Starting values fixed: gamma =  19 , th =  8.22092 ; SSE =  112.7901
Optimization algorithm converged
Optimized values fixed for regime 2  : gamma =  18.99999 , th =  8.246298
Testing for addition of regime 3.
  Estimating gradient matrix...
  Done. Computing the test statistic...
  Done. Regime 3 is needed (p-Value =  0.04979377 ).
Adding regime  3 .
Fixing good starting values for regime  3 ...
Reordering regimes...
Estimating parameters of regime 3 ...
Optimized values fixed for regime  3 : gamma =  39.99453 , th =  8.309892
*** Convergence problem. Code:  10
Optimized linear values:
0.3213105 1.298376 -0.3471048 -0.04305855
-27711956 2790671 946151.6 -269514.3
27711962 -2790672 -946151 269513.9
Ok.
Testing for addition of regime  4 .
  Estimating gradient matrix...
  Computing the test statistic...
Regime  4  is needed (p-Value =  0.03521662 ).
Finished building a MRSTAR with  3  regimes
> mod2.star
Non linear autoregressive model
Multiple regime STAR model
Regime  1 :
    Linear parameters: 0.3213105, 1.2983759, -0.3471048, -0.0430585
Regime  2 :
    Linear parameters: -27711956.1639132, 2790671.1486042, 946151.575765, -269514.305401
    Non-linear parameters:
19.0446151, 8.3207742
Regime  3 :
    Linear parameters: 27711962.1600718, -2790671.5917488, -946151.0472262, 269513.9229495
    Non-linear parameters:
39.9945317, 8.3098917

前辈您好，下面的程序是我对估计出的gamma做bootstrap，麻烦您看下有什么问题没有，附件里是我的说明，数据仍为quartc，还有又没有什么其他简单的包实现我想要的结果吗？谢谢您了!
library(tsDyn)
svpdx =read.table("quartc.txt", header = TRUE);
x =svpdx$qcpi
mod =star(x, m=5,d=1,thDelay=1,noRegimes=2,control=list(maxit=3000))
re=mod$residuals
i=seq(1,88,1)
n= length(i)
B=9999   
for  (b in 1:B) {
rr=c(6:n)
for (t in 6 :n){
number=ceiling(83*runif(1))
rr[t]=(re[number]-mean(re))/ sqrt(83/78)}
y=c(x[1:5], 6 :n)
for  (t in 6:n)
{y[t]=0.18831613+1.43411762*y[t-1]-0.31978565*y[t-2]-0.13549745*y[t-3]-0.19130543*y[t-4]+0.16428765*y[t-5]+(14.42516390-0.79064384*y[t-1]+0.70401255*y[t-2]+0.01523312*y[t-3]-0.57587624*y[t-4]-
0.08938476*y[t-5])/(1+exp(-23.99999222*(y[t-2]-13.99886588)))+rr[t]}
z=y [1:n]
rs=star(z, m=5,d=1,thDelay=1,noRegimes=2,control=list(maxit=3000))
gamma=rs$coefficients[13]}
qq=quantile(gamma,probs=seq(0,1,0.05))

2# epoh

#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=2,control=list(maxit=3000))

phi1=mod$model.specific$phi1
phi1
#>         [,1]            [,2]              [,3]               [,4]               [,5]               [,6]
#[1,]  0.1883161    1.4341176  -0.3197857 -0.13549745 -0.1913054  0.16428765
#[2,] 14.4251639 -0.7906438   0.7040126  0.01523312  -0.5758762  -0.08938476
phi2=mod$model.specific$phi2
phi2
#>        [,1]       [,2]
#[1,] 23.99999 13.99887
re=mod$residuals
n=length(x)    #88
n1=length(re)  #83
number=0*seq(1:n1)
B=300   
gamma=0*seq(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=re[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],phi2[2]) + rr[t-5]}
  
#{y[t]=0.18831613+1.43411762*y[t-1]-0.31978565*y[t-2]-0.13549745*y[t-3]-0.19130543*y[t-4]+0.16428765*y[t-5]+
#(14.42516390-0.79064384*y[t-1]+0.70401255*y[t-2]+0.01523312*y[t-3]-0.57587624*y[t-4]-
#0.08938476*y[t-5])/(1+exp(-23.99999222*(y[t-2]-13.99886588)))+re[t-5]}  #generate new y
z=y
rs=star(z, m=5,d=1,thDelay=1,noRegimes=2,trace=FALSE,control=list(maxit=3000))
if (is.null(rs$coefficients[13]) == TRUE) gamma=NA  else
      gamma=rs$coefficients[13]
}
gamma=gamma[!is.na(gamma)]          # remove NAs
length(gamma)     # 66
gamma
> gamma
[1] 40.001540 15.999911 14.970791  9.952839 40.000530  9.997871 40.001006
[8] 20.001059 14.038017  9.989034 11.900955 39.999946  9.999880 40.000001
[15] 16.999996  9.999980 40.056063 14.000000 11.000092 69.168903 15.000002
[22] 13.987616  9.699027 10.001973 40.000462  9.997138 28.000000 40.000186
[29] 40.008200 17.000007 19.001154  9.060692 13.999967 40.000230 18.000005
[36]  9.972111  9.999646 40.000281 31.000172 14.986656  9.997438 10.043781
[43] 20.000006 40.000001 17.000004  9.998270 40.000228 12.107126 40.000484
[50] 32.000004 21.001668  1.349424  9.997498  3.316196 40.000015 31.999999
[57] 40.000263 27.000000 23.000562 12.999567  9.964640 15.316762 40.270511
[64] 40.000000 33.000291 11.000011
qq=quantile(gamma,probs=seq(0,1,0.05))
qq
> qq
       0%        5%       10%       15%       20%       25%       30%
1.349424  9.762480  9.980572  9.997483  9.999646 10.012425 11.450523
      35%       40%       45%       50%       55%       60%       65%
13.740604 14.038017 15.079192 17.000000 18.750867 21.001668 28.750043
      70%       75%       80%       85%       90%       95%      100%
32.500148 40.000001 40.000186 40.000268 40.000507 40.006535 69.168903

前辈您好!我的意思你肯定已经明白了（我想构造LSTAR模型估计的所有系数的95%的置信区间，循环估计1000次），下面的程序可以实现，但是所需时间很长且好多都通不过线性检验(按说不应该出现这样的情况啊，由非线性方程产生的数据应该具有非线性性啊，即使是出现变异也不应该会有这样大的变异啊，开始时我怀疑过我的程序，但是我看不出程序那里有问题)，想问下您如果存在通不过检验的我构造置信区间的话还有意义吗？1000次的过程中还要不要估计那些通不过线性检验的模型呢？谢谢您了！
library(tsDyn)
svpdx =read.table("quartc.txt", header = TRUE);
x =svpdx$qcpi
mod =star(x, m=5,d=1,thDelay=1,noRegimes=2,control=list(maxit=3000))
re=mod$residuals

i=seq(1,88,1)
n= length(i)
B=1000
gamma=c(1:B)
coe1=c(1:B)
coe2=c(1:B)
coe3=c(1:B)
coe4=c(1:B)
coe5=c(1:B)
coe6=c(1:B)
coe7=c(1:B)
coe8=c(1:B)
coe9=c(1:B)
coe10=c(1:B)
coe11=c(1:B)
coe12=c(1:B)
coe14=c(1:B)
for(b in 1:B) {
rr=c(6:n)
for (t in 6:n){
number=ceiling(83*runif(1))
rr[t]=(re[number]-mean(re))*sqrt(83/78)}
y=c(x[1:5], 6 :n)
for
(t in 6:n)
{y[t]=0.18831613+1.43411762*y[t-1]-0.31978565*y[t-2]-0.13549745*y[t-3]-0.19130543*y[t-4]+0.16428765*y[t-5]+(14.42516390-0.79064384*y[t-1]+0.70401255*y[t-2]+0.01523312*y[t-3]-0.57587624*y[t-4]-0.08938476*y[t-5])/(1+exp(-23.99999222*(y[t-2]-13.99886588)))+rr[t]}
z=y [1:n]
rs=star(z, m=5,d=1,thDelay=1,noRegimes=2,sig=1,control=list(maxit=3000))
gamma=rs$coefficients[13]
coe1=rs$coefficients[1]
coe2=rs$coefficients[2]
coe3=rs$coefficients[3]
coe4=rs$coefficients[4]
coe5=rs$coefficients[5]
coe6=rs$coefficients[6]
coe7=rs$coefficients[7]
coe8=rs$coefficients[8]
coe9=rs$coefficients[9]
coe10=rs$coefficients[10]
coe11=rs$coefficients[11]
coe12=rs$coefficients[12]
coe14=rs$coefficients[14]
}

gg=gamma[1:B]
co1= coe1[1:B]
co2= coe2[1:B]
co3= coe3[1:B]
co4= coe4[1:B]
co5= coe5[1:B]
co6= coe6[1:B]
co7= coe7[1:B]
co8= coe8[1:B]
co9= coe9[1:B]
co10= coe10[1:B]
co11= coe11[1:B]
co12= coe12[1:B]
co14= coe14[1:B]

qq=quantile(gg,probs=seq(0,1,0.05))
cco1= quantile(co1,probs=seq(0,1,0.05))
cco2= quantile(co2,probs=seq(0,1,0.05))
cco3= quantile(co3,probs=seq(0,1,0.05))
cco4= quantile(co4,probs=seq(0,1,0.05))
cco5= quantile(co5,probs=seq(0,1,0.05))
cco6= quantile(co6,probs=seq(0,1,0.05))
cco7= quantile(co7,probs=seq(0,1,0.05))
cco8= quantile(co8,probs=seq(0,1,0.05))
cco9= quantile(co9,probs=seq(0,1,0.05))
cco10= quantile(co10,probs=seq(0,1,0.05))
cco11= quantile(co11,probs=seq(0,1,0.05))
cco12= quantile(co12,probs=seq(0,1,0.05))
cco14= quantile(co14,probs=seq(0,1,0.05))
cco1
cco2
cco3
cco4
cco5
cco6
cco7
cco8
cco9
cco10
cco11
cco12
cco14

63# epoh

前辈您好！我做bootstrap构造系数的置信区间时，我看到其他人都将残差标准化，如果不标准化的话做出来的结果会不会不可靠呢？还有由非线性方程产生的数据应该具有非线性性啊，即使是出现变异也不应该会有这样大的变异啊，如果剔除通不过检验的我构造置信区间的话还有意义吗？谢谢您了啊！ 64# 南冰

Bootstrap confidence intervals
请参考;
   http://www.colorado.edu/geograph ... b3/html/lab3.html#4

由于这个模型比较特殊,
会先做Linearity testing
所以会剔除条件不符者,
留下的都符合条件.
我个人认为还是具有代表性.

只是受限于原代码
maxGamma <- 40;
minGamma <- 10;
rateGamma <- 5;
范围就在10 ~40 之间,
要数据漂亮,可能要放宽范围,
0 ~100之间,但速度会变慢!

前辈您好!为什么我运行你的程序出问题了呢？麻烦您帮我看下，谢谢您了！
> library(tsDyn)
载入需要的程辑包：mgcv
This is mgcv 1.7-3. For overview type 'help("mgcv-package")'.
载入需要的程辑包：Matrix
载入需要的程辑包：lattice
载入程辑包：'Matrix'
The following object(s) are masked from 'package:base':
    det
载入需要的程辑包：snow
载入程辑包：'snow'
The following object(s) are masked from 'package:base':
    enquote
载入需要的程辑包：mnormt
载入需要的程辑包：foreach
载入需要的程辑包：iterators
载入需要的程辑包：codetools
foreach: simple, scalable parallel programming from REvolution Computing
Use REvolution R for scalability, fault tolerance and more.
http://www.revolution-computing.com
载入需要的程辑包：MASS
> svpdx =read.table("quartc.txt", header = TRUE);
> x =svpdx$qcpi
> mod =star(x, m=5,d=1,thDelay=1,noRegimes=2,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
Finished building a MRSTAR with 2 regimes
>
> phi1=mod$model.specific$phi1
> phi1
           [,1]       [,2]       [,3]        [,4]       [,5]        [,6]
[1,]  0.1883161  1.4341176 -0.3197857 -0.13549745 -0.1913054  0.16428765
[2,] 14.4251639 -0.7906438  0.7040126  0.01523312 -0.5758762 -0.08938476
> phi2=mod$model.specific$phi2
> phi2
         [,1]     [,2]
[1,] 23.99999 13.99887
> re=mod$residuals
> n=length(x)    #88
> n1=length(re)  #83
> number=0*seq(1:n1)
> B=300   
> gamma=0*seq(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=re[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],phi2[2]) + rr[t-5]}
+   
+ #{y[t]=0.18831613+1.43411762*y[t-1]-0.31978565*y[t-2]-0.13549745*y[t-3]-0.19130543*y[t-4]+0.16428765*y[t-5]+
+ #(14.42516390-0.79064384*y[t-1]+0.70401255*y[t-2]+0.01523312*y[t-3]-0.57587624*y[t-4]-
+ #0.08938476*y[t-5])/(1+exp(-23.99999222*(y[t-2]-13.99886588)))+re[t-5]}  #generate new y
+ z=y
+ rs=star(z, m=5,d=1,thDelay=1,noRegimes=2,trace=FALSE,control=list(maxit=3000))
+ if (is.null(rs$coefficients[13]) == TRUE) gamma=NA  else
+       gamma=rs$coefficients[13]
+ }
lstar: missing value during computations
lstar: missing value during computations
lstar: missing value during computations
lstar: missing value during computations
lstar: missing value during computations
警告信息：
1: In plogis(q, location, scale, lower.tail, log.p) : 产生了NaNs
2: In plogis(q, location, scale, lower.tail, log.p) : 产生了NaNs
3: In plogis(q, location, scale, lower.tail, log.p) : 产生了NaNs
4: In plogis(q, location, scale, lower.tail, log.p) : 产生了NaNs
5: In plogis(q, location, scale, lower.tail, log.p) : 产生了NaNs
> gamma=gamma[!is.na(gamma)]          # remove NAs
> length(gamma)     # 66
[1] 0
> length(gamma)     # 66
[1] 0
> qq=quantile(gamma,probs=seq(0,1,0.05))
> qq
  0%   5%  10%  15%  20%  25%  30%  35%  40%  45%  50%  55%  60%  65%  70%  75%  80%  85%  90%  95% 100%
  NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA



63# epoh

b先设小,b=50
多试几次

前辈您好，再请教你个问题啊！我在R的官网上下载到的包有两种，一种是包含源程序的，另一种是不包含源程序的，我选择本地安装时不包含源程序的package能安装成功，但是含有有原程序的package安装时老是提示错误，请问我如何安装自己改过原始CODE的package呢？谢谢您了！
R version 2.11.1 (2010-05-31)
Copyright (C) 2010 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R是自由软件，不带任何担保。
在某些条件下你可以将其自由散布。
用'license()'或'licence()'来看散布的详细条件。
R是个合作计划，有许多人为之做出了贡献.
用'contributors()'来看合作者的详细情况
用'citation()'会告诉你如何在出版物中正确地引用R或R程序包。
用'demo()'来看一些示范程序，用'help()'来阅读在线帮助文件，或
用'help.start()'通过HTML浏览器来看帮助文件。
用'q()'退出R.
[原来保存的工作空间已还原]
> utils:::menuInstallLocal()
错误于gzfile(file, "r") : 无法打开链结
此外: 警告信息：
1: In unzip(zipname, exdir = dest) : 从zip文件中抽取1时出了错
2: In gzfile(file, "r") :
  无法打开压缩文件'tsDyn_0.7-40.tar.gz/DESCRIPTION'，可能是因为'No such file or directory'
>

68# epoh

source code 是当你研究时,发现package
需要修改代码,以符合你的各项需求时,
供你修改用的.

修改完的source code需要经compile
re-build package成xxx.zip(ex:tsDyn_0.7-40.zip )
才能由本地安装

build package我已在58楼提过.

build package 需要安装

1.Rtools*.exe
        http://www.murdoch-sutherland.com/Rtools/index.html

主要含

*The command line tools (in Rtools*.exe)

*The MinGW toolchain to compile C, Fortran and C++.

2.LaTeX
        http://www.miktex.org/

你可选择basic version也可选择 complete version

我的电脑,安装完之后,路径如下:

PATH=c:\Rtools\bin;c:\Rtools\MinGW\bin;c:\MiKTeX\miktex\bin;

C:\Program Files\R\R-2.12.2\bin\i386;
请注意:你的会跟我不同,尤其是R

R的部分需要自行加入

详细步骤请参考:

Making R Packages Under Windows.pdf
  https://bbs.pinggu.org/thread-920948-1-1.html