楼主: 飘零的枫叶
51966 206

[问答] 用lingo软件如何对DEA-DA模型进行编程 [推广有奖]

121
ywh19860616 发表于 2011-7-30 09:27:22
用load命令,如数据名称为data
load data.dat

那么在左边Workspace会显示数据,打开就能看了
已有 2 人评分学术水平 热心指数 信用等级 收起 理由
zhangtao + 5 + 5 + 5 热心帮助其他会员
epoh + 1 + 1 + 1 我很赞同

总评分: 学术水平 + 6  热心指数 + 6  信用等级 + 6   查看全部评分

一份耕耘,一份收获。

122
epoh 发表于 2011-8-1 09:58:20
打开sp500.DAT文件
  load sp500.DAT

打开data.mat文件
  load data.mat

%%%%%%%%%%data_analysis.m
data.mat 有4个 variables
% p: average labor productivity
% u: unemployment rate
% v: vacancies posted
% y: gdp
Ch2007.xls也有4个 variables
可以套用,但未必合适!
因为作者是依其需要取对数
并增加一变量theta=v./u;

请自行斟酌.

%%%%%%%
X=xlsread('Ch2007.xls','B2:E57');
p=X(:,1);
u=X(:,2);
y=X(:,3);
v=X(:,4);

N=56;
smooth=1600; %  hpfilter smoothing papameter
theta=v./u;

logp=log(p);
logu=log(u);
logy=log(y);
logv=log(v);
logtheta=log(theta);

hp_p=hpfilter(logp,smooth);
hp_u=hpfilter(logu,smooth);
hp_y=hpfilter(logy,smooth);
hp_v=hpfilter(logv,smooth);
hp_theta=hpfilter(logtheta,smooth);

pdev=logp-hp_p;
udev=logu-hp_u;
ydev=logy-hp_y;
vdev=logv-hp_v;
thetadev=logtheta-hp_theta;

stddev=std([pdev, udev, vdev, thetadev, ydev]);
corrall=corrcoef([udev, vdev, thetadev]);

num=corrcoef([pdev(1:N-1), pdev(2:N)]);
ap=num(1,2);
num=corrcoef([udev(1:N-1), udev(2:N)]);
au=num(1,2);
num=corrcoef([vdev(1:N-1), vdev(2:N)]);
av=num(1,2);
num=corrcoef([thetadev(1:N-1), thetadev(2:N)]);
atheta=num(1,2);
num=corrcoef([ydev(1:N-1), ydev(2:N)]);
ay=num(1,2);

num=corrcoef([pdev(2:N),udev(1:N-1)]);
corrpu(1,1)=num(1,2);
num=corrcoef([pdev(1:N), udev(1:N)]);
corrpu(1,2)=num(1,2);
num=corrcoef([pdev(1:N-1), udev(2:N)]);
corrpu(1,3)=num(1,2);
num=corrcoef([pdev(1:N-2), udev(3:N)]);
corrpu(1,4)=num(1,2);
num=corrcoef([pdev(1:N-3), udev(4:N)]);
corrpu(1,5)=num(1,2);

num=corrcoef([pdev(2:N), vdev(1:N-1)]);
corrpv(1,1)=num(1,2);
num=corrcoef([pdev(1:N), vdev(1:N)]);
corrpv(1,2)=num(1,2);
num=corrcoef([pdev(1:N-1), vdev(2:N)]);
corrpv(1,3)=num(1,2);
num=corrcoef([pdev(1:N-2), vdev(3:N)]);
corrpv(1,4)=num(1,2);
num=corrcoef([pdev(1:N-3), vdev(4:N)]);
corrpv(1,5)=num(1,2);

num=corrcoef([pdev(2:N), thetadev(1:N-1)]);
corrptheta(1,1)=num(1,2);
num=corrcoef([pdev(1:N), thetadev(1:N)]);
corrptheta(1,2)=num(1,2);
num=corrcoef([pdev(1:N-1), thetadev(2:N)]);
corrptheta(1,3)=num(1,2);
num=corrcoef([pdev(1:N-2), thetadev(3:N)]);
corrptheta(1,4)=num(1,2);
num=corrcoef([pdev(1:N-3), thetadev(4:N)]);
corrptheta(1,5)=num(1,2);

num=corrcoef([ydev(2:N), pdev(1:N-1)]);
corryp(1,1)=num(1,2);
num=corrcoef([ydev(1:N), pdev(1:N)]);
corryp(1,2)=num(1,2);
num=corrcoef([ydev(1:N-1), pdev(2:N)]);
corryp(1,3)=num(1,2);

num=corrcoef([ydev(2:N), udev(1:N-1)]);
corryu(1,1)=num(1,2);
num=corrcoef([ydev(1:N), udev(1:N)]);
corryu(1,2)=num(1,2);
num=corrcoef([ydev(1:N-1), udev(2:N)]);
corryu(1,3)=num(1,2);

num=corrcoef([ydev(2:N), vdev(1:N-1)]);
corryv(1,1)=num(1,2);
num=corrcoef([ydev(1:N), vdev(1:N)]);
corryv(1,2)=num(1,2);
num=corrcoef([ydev(1:N-1), vdev(2:N)]);
corryv(1,3)=num(1,2);

num=corrcoef([ydev(2:N), thetadev(1:N-1)]);
corrytheta(1,1)=num(1,2);
num=corrcoef([ydev(1:N), thetadev(1:N)]);
corrytheta(1,2)=num(1,2);
num=corrcoef([ydev(1:N-1), thetadev(2:N)]);
corrytheta(1,3)=num(1,2);

已有 1 人评分学术水平 热心指数 信用等级 收起 理由
zhangtao + 5 + 5 + 5 好的意见建议

总评分: 学术水平 + 5  热心指数 + 5  信用等级 + 5   查看全部评分

123
zhangtao 发表于 2011-8-1 13:59:17
epoh 发表于 2011-8-1 09:58
打开sp500.DAT文件
  load sp500.DAT
epoh老师,附件中的eviews程序,我估计出的结果有一部分是:NA,您看一下是什么原因?
为什么?如何修改?
非常感谢!
LogL: LL_L   
Method: Maximum Likelihood (Marquardt)   
Date: 08/01/11   Time: 13:47   
Sample: 3 195   
Included observations: 193   
Evaluation order: By observation   
Estimation settings: tol= 1.0e-05, derivs=accurate numeric   
Initial Values: MU(1)=0.13502, OMEGA(1)=1.45022, ALPHA(1)=0.04387,   
        BETA(1)=0.70089, TDF(1)=3.00000   
Convergence achieved after 3 iterations   
WARNING: Singular covariance - coefficients are not unique   
   
Coefficient Std. Error z-Statistic Prob.  
   
MU(1) 0.212837 NA NA NA
OMEGA(1) -3.361117 NA NA NA
ALPHA(1) -0.034662 NA NA NA
BETA(1) 1.556563 NA NA NA
TDF(1) 7.206290 NA NA NA
   
Log likelihood -436.6545     Akaike info criterion  4.576730
Avg. log likelihood -2.262458     Schwarz criterion  4.661256
Number of Coefs. 5     Hannan-Quinn criter.  4.610961
   


数学好就是要天天学

124
zhangtao 发表于 2011-8-1 14:33:08
x=c(5,4,9,6,21,17,11,20,7,10,21,15,13,16,8)
n=length(x)
boot=10000
s2=NA*seq(1:boot)
for(i in 1:boot)
{
xsamp=sample(x,n,replace = T)
s2=sum((xsamp-mean(xsamp))^2)/(n-1)
}
var(s2)
#[1] 57.49942
本文来自: 人大经济论坛 S-Plus&R专版 版,详细出处参考: https://bbs.pinggu.org/forum.php? ... 1&from^^uid=11232

2# epoh

s2=NA*seq(1:boot)     epoh老师,您好!这一行是什么意思?NA并没有赋值,NA乘以seq序列是做什么用的?
非常感谢!

本文来自: 人大经济论坛 S-Plus&R专版 版,详细出处参考: https://bbs.pinggu.org/forum.php? ... 1&from^^uid=11232
数学好就是要天天学

125
epoh 发表于 2011-8-1 18:30:46

一个变量使用前要先定义

不管是

scale      s=5;

vector     v=seq(1:6)

matrix     m=matrix(v,3,2);

boot=100

我定义    s2=NA*seq(1:boot)

主要是方便除错

如果还有NA,就表示程序有误.

e-views code

   e-views code.rar (4.16 KB)

已有 1 人评分学术水平 热心指数 信用等级 收起 理由
zhangtao + 5 + 5 + 5 对epoh大师致以最真诚的感谢!!!

总评分: 学术水平 + 5  热心指数 + 5  信用等级 + 5   查看全部评分

126
zhangtao 发表于 2011-8-1 22:02:05
epoh 发表于 2011-8-1 18:30
一个变量使用前要先定义不管是scale      s=5;vector     v=seq(1:6)matrix     m=matrix(v,3,2); boot=100 ...
s2=sum((xsamp-mean(xsamp))^2)/(n-1)

本文来自: 人大经济论坛 Matlab及其他计量软件专版 版,详细出处参考: https://bbs.pinggu.org/forum.php?mod=viewthread&tid=1127506&page=13&from^^uid=11232

epoh老师,您好!
那s2定义了两次,您能解释一下吗?我还是不明白。
s2=NA*seq(1:boot) 这一行是如何纠错的?
本文来自: 人大经济论坛 Matlab及其他计量软件专版 版,详细出处参考: https://bbs.pinggu.org/forum.php?mod=viewthread&tid=1127506&page=13&from^^uid=11232

另外, zhang.rar (173.33 KB) 附件中的pftest_d.m运行时总是提示以下错误:
??? Error using ==> var
W must be a vector of nonnegative weights, or a scalar 0 or 1.
Error in ==> pftest_d at 17
results = var(y,nlag);
您看是什么问题?如何解决?
另外,epoh老师,我的matlab7没有var.m这个文件,您能把您的给我传一份吗?
非常感谢!


>> randn('state',0)                % Start from a known state.
x = randn(1000,1);              % 1000 Gaussian deviates ~ N(0,1).
y = filter(1,[1 -0.6 0.08],x);  % Create a stationary AR(2)
                                % process.
[PartialACF,Lags,Bounds] = parcorr(y,[],2); % Compute the
                                % partial ACF with 95 percent
                                % confidence.
[Lags,PartialACF]
Bounds
parcorr(y,[],2)            % Use the same example, but plot
                           % the partial ACF sequence with
                           % confidence bounds.
??? Error using ==> lagmatrix
lagmatrix: wrong # of input arguments
Error in ==> parcorr at 151
X          =  lagmatrix(Series , [1:nLags]);
>>
这段代码出错,为什么?如何修改?
非常感谢!

数学好就是要天天学

127
zhangtao 发表于 2011-8-1 22:15:43
epoh老师,您好!
还有这个ecm_d.m,运行时提示:
??? Error using ==> fopen
The file mode for fopen must contain exactly one of the modes 'r', 'w' or 'a'.

Error in ==> ecm_d at 33
fid = fopen('ecm.out','wr');

>>
这又如何解决?
非常感谢!
数学好就是要天天学

128
epoh 发表于 2011-8-2 09:41:12
jplv7 toolbox有一些小地方没更正,
恰好提供自我训练的机会

1.Error using ==> var
  请键入help var
  你会发现用到的var.m是
  matlab 自带的var(),是计算Variance.
  想要performs vector autogressive estimation
  作者已由var(),改为vare()
  所以results = var(y,nlag);
  请改results = vare(y,nlag);

2.同理.Error using ==> lagmatrix
  请键入help lagmatrix
  你会发现用到的lagmatrix.m是
  jplv7\Ucsd_garch\Garch\lagmatrix.m
  非matlab 自带的lagmatrix()
  所以请先将jplv7\Ucsd_garch\Garch\lagmatrix.m
  改个名字

3.Error using ==> fopen
  The file mode for fopen must contain exactly
  one of the modes 'r', 'w' or 'a'.
  请将fid = fopen('ecm.out','wr');
  改为fid = fopen('ecm.out','w');


#######
我用简单的bootstrap=2说明
第二次用到s2
是把计算出来的值依序放进s2[1],s2[2],..
然后计算var(s2)
x=c(5,4,9,6,21,17,11,20,7,10,21,15,13,16,8)
n=length(x)
boot=2
s2=NA*seq(1:boot)
for(i in 1:boot)
{
xsamp=sample(x,n,replace = T)
s2[i]=sum((xsamp-mean(xsamp))^2)/(n-1)
}
var(s2)

[1] "xsamp=" 17 16  5  6  5  9 13 15 11 20  6 21  5  4 21
[1] "s2="   40.54286       NA

[1] "xsamp="  5  5  5  7 21 21  5  5  8 21 17  7  7 15  9
[1] "s2="   40.54286 42.12381
已有 1 人评分学术水平 热心指数 信用等级 收起 理由
ywh19860616 + 5 + 5 + 5 好的建议

总评分: 学术水平 + 5  热心指数 + 5  信用等级 + 5   查看全部评分

129
ywh19860616 发表于 2011-8-2 09:58:26
可以打开jplv7目录下var_bar,然后点击contents.html,这个网页介绍了各个函数的用法

只是现在没有panel data的VAR程序
一份耕耘,一份收获。

130
zhangtao 发表于 2011-8-2 10:16:25
> x=c(5,4,9,6,21,17,11,20,7,10,21,15,13,16,8)
> n=length(x)
> boot=2
> s2=NA*seq(1:boot)
> for(i in 1:boot)
+ {
+ xsamp=sample(x,n,replace = T)
+ s2=sum((xsamp-mean(xsamp))^2)/(n-1)
+ }
> var(s2)
[1] 30.94222
>
> 本文来自: 人大经济论坛 Matlab及其他计量软件专版 版,详细出处参考: https://bbs.pinggu.org/forum.php? ... 3&from^^uid=11232
错误: 意外的符号在"本文来自: 人大经济论坛 Matlab及其他计量软件专版"里
> xsample
错误: 找不到对象'xsample'
> xsamp
[1]  8 21 15 21  9 10 21 21  5 11 21 10  9 15 15
> n
[1] 15
> s2
[1] 24.68571 32.55238
> s2$1
错误: 意外的数值量在"s2$1"里
> s2[1]
[1] 24.68571
> s2[2]
[1] 32.55238
> s2[3]
[1] NA
>
> s2[0]
numeric(0)
> s2=NA*seq(1:boot)
> s2
[1] NA NA
> s2=sum((xsamp-mean(xsamp))^2)/(n-1)
> s2
[1]       NA 32.55238
> s2
[1]       NA 32.55238
> s2[2]
[1] 32.55238
> s2[1]
[1] NA
> mean(xsamp
+
+ mean(xsamp)
错误: 意外的符号 于
"
mean"
> xsamp
[1]  8 21 15 21  9 10 21 21  5 11 21 10  9 15 15
> mean(xsamp)
[1] 14.13333
> x=c(5,4,9,6,21,17,11,20,7,10,21,15,13,16,8)
> n=length(x)
> boot=2
> s2=NA*seq(1:boot)
> for(i in 1:boot)
+ {
+ xsamp=sample(x,n,replace = T)
+ s2=sum((xsamp-mean(xsamp))^2)/(n-1)
+ }
> var(s2)
[1] 10.44898
> s2
[1] 39.83810 44.40952
epoh老师,您好!
    我还是不明白:> s2=NA*seq(1:boot)
> s2
[1] NA NA

> n
[1] 15


> s2=sum((xsamp-mean(xsamp))^2)/(n-1)
> s2
[1]       NA 32.55238
> s2
[1]       NA 32.55238

第二次运行程序:

> x=c(5,4,9,6,21,17,11,20,7,10,21,15,13,16,8)
> n=length(x)
> boot=2
> s2=NA*seq(1:boot)
> for(i in 1:boot)
+ {
+ xsamp=sample(x,n,replace = T)
+ s2=sum((xsamp-mean(xsamp))^2)/(n-1)
+ }
> var(s2)
[1] 10.44898
> s2
[1] 39.83810 44.40952
第一次运行结果和第二次运行结果相矛盾,不一致,
另外,这一行: s2=NA*seq(1:boot)
您说是为了纠错,那纠错的内在机制是什么?
我还是不明白。为什么第一次和第二次运行的结果会矛盾?
非常感谢!

数学好就是要天天学

您需要登录后才可以回帖 登录 | 我要注册

本版微信群
加好友,备注cda
拉您进交流群
GMT+8, 2026-2-3 07:38