全局变量的使用问题

5588

收藏 2011-10-10

首先在此感谢epoh的热心帮助。我想使用全局变量，是因为把函数结构传递到这个变量上，以后不需要每次都进行选择判断。代码1和代码2唯一的区别是：代码1，我每个子函数都定义了TR，代码2，我只在主函数定义全局变量TR，子函数申明全局变量TR。我检查了子函数TR的传递情况，确实是能够传递过去，但不知道为什么结果相差那么大？代码1结果比较接近真实值。另外一个问题是，初始值的选择问题，我用的ols回归和格点搜索进行选择，但monte caro模拟发现这个方法很失败？不知道大家是怎么选择初始值的
我把代码和测试都传上来，希望大家多给予指导指导！

附件列表

code.txt

大小:7.88 KB

马上下载

扫码加我拉你入群

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

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

全部回复

hongri200890

2011-10-10 20:22:33

怎么到经济论坛上问计算机啊……

扫码加我拉你入群

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

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

epoh

2011-10-10 21:38:14

老兄

1.这个程序,有些乱吧

乱到不知从何改起

TR是一个 Inline function

你为甚么不分别定义

LSTAR & ESTAR的 Transition function

要用时,带入参数呼叫即可

就如R底下的定义

#Transition function

#y: variable

#g: smoothing parameter

#c: threshold value

G <- function(y, g, th)

plogis(y, th, 1/g)

你先自行修改应能解决

你所谓全局的问题

2.c: threshold value要定范围

否则会产生NA,就如R底下

# 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;

3.如果你是做这个专题

请考虑用R package "tsDyn"

虽然目前只有lstar

但他采grid search for starting values

精准许多

扫码加我拉你入群

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

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

kerrydu

2011-10-10 21:57:46

非常感谢epoh兄，为了管理方便我把很多子程序都写到一起，看起来确实比较乱。没有分别定义ESTAR&LSTAR是因为后面每次碰到都得用if进行判断，我想一开始就判断，然后把状态转移函数传递到一个TR，以后就不要再做选择判断，这样效率会高很多，特别是优化的迭代的时候。我想不明白的是这两个代码几乎一样，结果相差那么大？唯一可能出现问题的是全局变量的传导，可我检查了子函数里TR的传递确实是成功的。
epoh兄是R高手呢，可惜我对R不太熟悉，以后还要向epoh多加请教呢

扫码加我拉你入群

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

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

kerrydu

2011-10-10 22:00:59

hongri200890 发表于 2011-10-10 20:22
怎么到经济论坛上问计算机啊……

呵呵，论坛上计算高手也很多哦

扫码加我拉你入群

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

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

hongri200890

2011-10-10 22:03:48

哦好吧……

扫码加我拉你入群

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

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

点击查看更多内容…

zhangtao

2011-10-11 09:35:31

epoh 发表于 2011-10-10 21:38
老兄1.这个程序,有些乱吧乱到不知从何改起 TR是一个 Inline function 你为甚么不分别定义 LSTAR & E ...

code.rar
大小:(9.32 KB)

马上下载

epoh老师，您好！
您看一下附件中程序运行后为什么会出现以下问题？如何修改？
非常感谢！
另外，建议楼主编程时多一些注释，这样更容易让其他人读懂。

运行楼主的的代码2
>> maincode
rpara =
0.3000
0.5000
-0.1000
0.7000
0.1200
0.5000
5.0000
0.8000
??? Error using ==> STAR_startval1
Too many input arguments.
Error in ==> STAR1 at 13
stval=STAR_startval1(para0,flag,y,X,q,0);%计算初始值，这个函数不好用，计算的初始值经常收敛到局部最小值
Error in ==> maincode at 3
para=STAR1(y,X,q,1,rpara)
运行楼主的代码1
>> maincode
rpara =
0.3000
0.5000
-0.1000
0.7000
0.1200
0.5000
5.0000
0.8000
??? Undefined command/function 'STAR_itransform'.
Error in ==> STARlik1 at 11
para=STAR_itransform(para);
Error in ==> fminunc at 182
f = feval(funfcn{3},x,varargin{:});
Error in ==> STAR1 at 18
[para,fv]=fminunc(@STARlik1,stval,options,flag,y,X,q,1);
Error in ==> maincode at 3
para=STAR1(y,X,q,1,rpara)
>>

扫码加我拉你入群

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

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

epoh

2011-10-11 11:23:20

建议你用fmincon

FMINCON finds a constrained minimum of a function of several variables

sig1,sig2,应该都是positive

c根据q,也要有一定范围

譬如:

stval =

0.3000

0.5000

-0.1000

0.7000

-2.1203

-0.6931

1.6094

-0.2231

%%%%

假设你的程序是对的

你要用global

请修改两个部分

1.程序一开始就要宣告

2.STAR2.m就不用宣告

%%%%%%
flag=1;
ep=randn(10000,1);data=zeros(10000,1);
data(1)=0.5;c=0.8;gamma=5;a1=[0.3;0.5];a2=[-0.1;0.7];sig1=0.12;sig2=0.5;
global TR;
if flag
TR=inline('1 ./ (1 + exp(-gamma*(q-c)))','gamma','q','c')
else
TR=inline('1- exp(-gamma*(q-c).^2)','gamma','q','c')
end

for k=2:10000
G1=TR(gamma,c,data(k-1));
data(k)=([1 data(k-1)]*a1+sig1*ep(k))*G1+([1 data(k-1)]*a2+sig2*ep(k))*(1-G1);
end

y=data(2:end);X=[ones(9999,1),data(1:end-1)];
q=data(1:end-1);% q----transiton varialbe
rpara=[a1;a2;sig1;sig2;gamma;c];% real parameter as one of start value
para1=STAR2(y,X,q,1,rpara)

fv =

-2.2161e+003

para1 =

0.2994
0.3076
0.2984
0.4207
0.4383
0.1249
5.4969
0.7390

扫码加我拉你入群

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

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

kerrydu

2011-10-11 11:32:00

我试过fmincon，问题是有时候最优值会落在你设定的边界上。sig1 sig2我已经通过exp变换转化成无约束问题了。STAR2就是主函数，在STAR2宣告应该没问题吧

扫码加我拉你入群

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

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

kerrydu

2011-10-11 12:12:50

zhangtao 发表于 2011-10-11 09:35
epoh老师，您好！
您看一下附件中程序运行后为什么会出现以下问题？如何修改？
非常感谢 ...

你是不是落下了子函数STAR——startval1了

扫码加我拉你入群

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

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

epoh

2011-10-11 16:08:30

不理解你的function STAR_tranform.m

为何gamma,c,也要取log,

一开始c starting value=-0.2231

很容易造成fv -1.3655e+003 +3.1413e+004i

stval =

0.3000
0.5000
-0.1000
0.7000
0.1200
0.5000
5.0000
0.8000

log(stval(end-3:end))

ans =

-2.1203
-0.6931
1.6094
-0.2231

fv =

-1.3655e+003 +3.1413e+004i

para1 =

0.6991
0.5624
0.3236
0.3302
-0.1698
-0.2111
2.6342
-3.6206

%%%%%%%%%%%

而在STARlik12.m

似乎没有转换回来

numxs=size(X,2)
gamma=para(end-1)
c=para(end)
sigma_1=para(end-3)
sigma_2=para(end-2)

phi1=para(1:numxs)
phi2=para(numxs+1:2*numxs)
mu1=X*phi1;
mu2=X*phi2;
G1=TR(gamma,q,c);
res=y-mu1.*G1-mu2.*(1-G1);

扫码加我拉你入群

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

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

zhangtao

2011-10-11 22:40:25

kerrydu 发表于 2011-10-11 12:12
你是不是落下了子函数STAR——startval1了

没有落下子函数，我用的是matlab7

扫码加我拉你入群

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

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

kerrydu

2011-10-13 11:39:50

epoh 发表于 2011-10-10 21:38
老兄1.这个程序,有些乱吧乱到不知从何改起 TR是一个 Inline function 你为甚么不分别定义 LSTAR & E ...

epoh兄，我这两天看了你给我推荐的R ，请教一下。
Usage
lstar(x,m,d=1,steps=d,series,mL,mH,thDelay,
th,gamma,trace=TRUE,control=list())

x timeseries
m,d,steps embeddingdimension,timedelay,forecastingsteps
series timeseriesname(optional)
mL autoregressiveorderfor’low’regime(dafult:m).Mustbe<=m
mH autoregressiveorderfor’high’regime(default:m).Mustbe<=m
这个m是什么意思
里面的例子，mod.lstar<-lstar(log10(lynx),m=2,mTh=c(0,1),control=list(maxit=3000))
mTh是一个向量，c（0,1）表示的是什么？我才接触R，不太明白

扫码加我拉你入群

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

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

epoh

2011-10-13 14:23:48

Original data:
yy=log10(lynx)---114 obs
[1] 2.429752 2.506505 2.767156  2.940018 3.168792 3.450403 3.594171 3.774006
.....
.....

[107] 3.186674 2.723456 2.685742 2.820858 3.000000 3.201397 3.424392 3.530968
#####
when m=2=mL=mH; d=1=steps;
model:
      yy    V1 V1.0 V1..1          V1.1             V1.0.1          V1..1.1
1 2.767156  1 2.506505 2.429752 4.477747e-05 1.122349e-04 1.087982e-04
2 2.940018  1 2.767156 2.506505 1.041599e-04 2.882267e-04 2.610773e-04
3 3.168792  1 2.940018 2.767156 1.828708e-03 5.376435e-03 5.060321e-03
4 3.450403  1 3.168792 2.940018 1.211853e-02 3.840111e-02 3.562871e-02
5 3.594171  1 3.450403 3.168792 1.318946e-01 4.550895e-01 4.179465e-01
.....
.....
108 2.820858  1 2.685742 2.723456 1.131569e-03 3.039102e-03 3.081778e-03
109 3.000000  1 2.820858 2.685742 7.476169e-04 2.108921e-03 2.007906e-03
110 3.201397  1 3.000000 2.820858 3.296504e-03 9.889513e-03 9.298970e-03
111 3.424392  1 3.201397 3.000000 2.317971e-02 7.420746e-02 6.953913e-02
112 3.530968  1 3.424392 3.201397 1.786306e-01 6.117009e-01 5.718673e-01
#########
mTh: coefficients for the lagged time series, to obtain the threshold variable
mTh: z[t] = x[t]mTh[1] + x[t-d]mTh[2] + .. + x[t-(m-1)d]mTh[m]
when mTh=c(0,1),z <- xx %*% mTh #threshold variable
xx=
      V1/0       V1/-1
  [1,] 2.506505 2.429752
  [2,] 2.767156 2.506505
  [3,] 2.940018 2.767156
  [4,] 3.168792 2.940018
  [5,] 3.450403 3.168792
  .....
  .....
  [108,] 2.685742 2.723456
  [109,] 2.820858 2.685742
  [110,] 3.000000 2.820858
  [111,] 3.201397 3.000000
  [112,] 3.424392 3.201397
thVar
[1] 2.429752 2.506505 2.767156 2.940018 3.168792 3.450403 3.594171 3.774006
.....
.....
[105] 3.553155 3.467608 3.186674 2.723456 2.685742 2.820858 3.000000 3.201397

扫码加我拉你入群

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

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

epoh

2011-10-14 08:56:14

做个总结

我把程序精简.

latsr script + lstarSSE.m (Objective function)

得出的结果跟 R package "tsDyn"接近.

para'
ans =  0.4883 1.2466 -0.3661 -1.0371 0.4237 -0.2518 11.0855 3.3396
####R package "tsDyn"
LSTAR model
Coefficients:
Low regime:
phi1.0    phi1.1    phi1.2
0.4870673  1.2465799 -0.3655721

High regime:
phi2.0 phi2.1 phi2.2
-1.0514786 0.4243230 -0.2488072

Smoothing parameter: gamma = 11

Threshold
Variable: Z(t) = + (0) X(t) + (1) X(t-1)

Value: 3.34
#############
%%%%%%%%%%%%%
%lstar script
%%%%%%%%%%%%%
load 'lynx.dat'
lynx=log10(lynx);
y=lynx(3:end,1);
X=[ones(length(y),1) lynx(2:(end-1),1) lynx(1:(end-2),1)];
q=lynx(1:(end-2),1);
%set initial values
gamma=20;
minTh=quantile(q,0.1) % 2.1353
maxTh=quantile(q,0.9) % 3.5788
c=(minTh+maxTh)/2
TR=inline('1 ./ (1 + exp(-gamma*(q-c)))','gamma','q','c')
XX=[X X(:,1).*TR(gamma,c,q) X(:,2).*TR(gamma,c,q) X(:,3).*TR(gamma,c,q)];
b=regress(y,XX)
init=[b; gamma; c]
lb=[-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,0,minTh ]
ub=[Inf,Inf,Inf,Inf,Inf,Inf,100,maxTh ]

[para]=fmincon(@(para)lstarSSE(para,y,X,q),init,[],[],[],[],lb,ub)

%%Objection function
%%%%%%%%%%%
%lstarSSE.m
%%%%%%%%%%%
function SSE=lstarSSE(para,y,X,q)
TR=inline('1 ./ (1 + exp(-gamma*(q-c)))','gamma','c','q');
numxs=size(X,2);
gamma=para(end-1);
c=para(end);
phi1=para(1:numxs);
phi2=para(numxs+1:2*numxs);
mu1=X*phi1;
mu2=X*phi2;
G1=TR(gamma,c,q);
res=y-mu1-mu2.*G1;
SSE=sum(res.^2);
%%%%%%%%%%%

lynx.dat

lynx.rar
大小:(386 Bytes)

马上下载

本附件包括：

lynx.dat

扫码加我拉你入群

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

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

kerrydu

2011-10-14 10:05:55

epoh 发表于 2011-10-14 08:56
做个总结我把程序精简.latsr script + lstarSSE.m (Objection function) 得出的结果跟 R package "tsDyn"接 ...

谢谢epoh兄

扫码加我拉你入群

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

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

扫码加我 拉你入群

扫码加我 拉你入群

扫码加我 拉你入群

扫码加我 拉你入群

扫码加我 拉你入群

扫码加我 拉你入群

扫码加我 拉你入群

扫码加我 拉你入群

扫码加我 拉你入群

扫码加我 拉你入群

扫码加我 拉你入群

扫码加我 拉你入群

扫码加我 拉你入群

扫码加我 拉你入群

扫码加我 拉你入群

扫码加我 拉你入群

分享

扫码加好友，拉您进群

扫码加我拉你入群

扫码加我拉你入群

扫码加我拉你入群

扫码加我拉你入群

扫码加我拉你入群

扫码加我拉你入群

扫码加我拉你入群

扫码加我拉你入群

扫码加我拉你入群

扫码加我拉你入群

扫码加我拉你入群

扫码加我拉你入群

扫码加我拉你入群

扫码加我拉你入群

扫码加我拉你入群

扫码加我拉你入群