全部版块 我的主页
论坛 数据科学与人工智能 数据分析与数据科学 R语言论坛
2011-3-29 19:46:21
Note:this modified package only for
         function star() m=4,ar(1),ar(3)
tsDyn_0.7-40.zip
大小:(1.32 MB)

 马上下载



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
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2011-3-29 21:19:14
前辈您好,下面的程序是我对估计出的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
附件列表

boostrap.doc

大小:39.5 KB

 马上下载

二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2011-3-30 20:52:39
#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
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2011-3-31 08:11:14
前辈您好!我的意思你肯定已经明白了(我想构造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
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

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

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2011-3-31 22:25:12
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之间,但速度会变慢!
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2011-3-31 22:29:51
前辈您好!为什么我运行你的程序出问题了呢?麻烦您帮我看下,谢谢您了!
> 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
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2011-3-31 22:57:27
b先设小,b=50
多试几次
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2011-4-1 09:36:35
前辈您好,再请教你个问题啊!我在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
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2011-4-1 12:15:50
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

二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2011-4-1 18:01:21
前辈您好!再请教你个问题啊。我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
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2011-4-1 21:49:46
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.
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2011-4-4 16:47:15
借个角落,来答覆,耽搁许久,
短信息询问系数标准差问题的朋友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
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2011-4-27 15:33:04
前辈,您好!我现在是想做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
附件列表

RSTAR.pdf

大小:1.18 MB

 马上下载

二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2011-4-28 13:39:16
##########
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

二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2011-5-10 13:43:39
75# epoh

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

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

祝福每一位热心人!

AZA AZA Fighting!
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2011-5-11 14:21:30
学习了 多谢楼主啊
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2011-5-11 16:32:08
epoh老师绝对是大牛啊,人品学术没的说啊,谢谢epoh老师啊! 76# O(∩_∩)O~!
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2011-5-11 20:27:41
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')

girf.bmp
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2011-5-16 22:29:20
前辈您好,非常感谢您给我编的程序,不过还得麻烦您,您写的程序有些地方我看不大懂,所以我不知道三区制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
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2011-5-17 19:19:09
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')
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2011-5-21 00:53:04
前辈您好!您可以把您在下面这个帖子的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
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2011-5-21 09:02:51
内含四个文件:
修改过 : lstar_new.R , isLinear.R
原始文件 : lstar.R , star.R


lstar files.rar
大小:(21.59 KB)

 马上下载

二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2011-5-21 13:57: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
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2011-5-21 15:10:03
你说的问题作者解决了
请下载最新版 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()可能忘记更新了,还是用我给你的.
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

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

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2011-5-21 20:10:38
estar.R
estar_new.rar
大小:(5.22 KB)

 马上下载

本附件包括:

  • estar.R

二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2011-5-21 21:41:32
前辈您好!麻烦您检查下发给我的这个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
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2011-5-22 10:11:54
重新检查一次
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
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2011-5-22 12:24:19
前辈您好!非常感谢您我一直以来对我的帮助,还有几个问题希望得到您的解答!
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
附件列表

m.txt

大小:382 Bytes

 马上下载

mzh.txt

大小:381 Bytes

 马上下载

二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

栏目导航
热门文章
推荐文章

分享

扫码加好友,拉您进群
各岗位、行业、专业交流群