用lingo软件如何对DEA-DA模型进行编程

刚开始学习DEA，时间比较紧迫需要尽快熟悉相关软件求解，首先谢过论坛朋友给予支持啊！

先请问，没有相关基础，什么软件最容易上手，lingo 还是 matlab ，这两个都不会呢，谢谢!

epoh 发表于 2011-8-6 14:08
data=xlsread('iotablecoe.xls','B2:K11');size(data) % 10 10           %你的数据有 10rows,10 columns%你 ...

20090191_data.rar (126.87 KB)

2011-8-7 09:20:57 上传

epoh老师，您好！
   为什么在运行附件中的数据和程序时会出现以下错误？如何解决？
非常感谢！
To get started, select MATLAB Help or Demos from the Help menu.
??? Error using ==> dlmread
The file 'C:\Users\Jason\Desktop\Medicare\HHsampnodum.csv' could not be opened because: No such file or directory
Error in ==> PartDHH4BS at 39
data = dlmread('C:\Users\Jason\Desktop\Medicare\HHsampnodum.csv');
>>

>  library(boot)#当然要先下载安装boot包
>  set.seed(200)
>  x=c(5,4,9,6,21,17,11,20,7,10,21,15,13,16,8)
>  x.boot <- boot(x,function(x,i)var(x[i]),R=999)
>  x.boot

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = x, statistic = function(x, i) var(x[i]), R = 999)


Bootstrap Statistics :
    original    bias    std. error
t1* 34.31429 -2.203737    7.322269
> a=x.boot$std.error**2
> var(a)
[1] NA
> a
numeric(0)
> a=x.boot$std.error^2
> a
numeric(0)
>
epoh老师，您好！
以上程序为什么a的值为numeric(0)？
a=x.boot$std.error^2这一行有什么错误？


> 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[i]=sum((xsamp-mean(xsamp))^2)/(n-1)
+ }
> var(s2)
[1] 58.19342
然后删除了s2=NA*seq(1:boot)
> x=c(5,4,9,6,21,17,11,20,7,10,21,15,13,16,8)
> n=length(x)
> boot=10000
>
> for(i in 1:boot)
+ {
+ xsamp=sample(x,n,replace = T)
+ s2[i]=sum((xsamp-mean(xsamp))^2)/(n-1)
+ }
错误于s2[i] = sum((xsamp - mean(xsamp))^2)/(n - 1) : 找不到对象's2'
> var(s2)
错误于inherits(x, "data.frame") : 找不到对象's2'
>
>

epoh老师，想问：s2=NA*seq(1:boot)这一行是不是对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)
> s2
[1] NA NA
> for(i in 1:boot)
+ {
+ xsamp=sample(x,n,replace = T)
+ s2[i]=sum((xsamp-mean(xsamp))^2)/(n-1)
+ }
> var(s2)
[1] 0.255102
> s2[3]
[1] NA
> s2[0]
numeric(0)
> s2
[1] 24.28571 25.00000
>
> 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)
> s2
[1] NA NA
> s2[1]
[1] NA
> for(i in 1:boot)
+ {
+ xsamp=sample(x,n,replace = T)
+ s2[i]=sum((xsamp-mean(xsamp))^2)/(n-1)
+ }
> var(s2)
[1] 123.0208
>
>

Error using ==> dlmread

  The file 'C:\Users\Jason\Desktop\Medicare\HHsampnodum.csv'

  could not be opened because: No such file or directory

  已经很明显告诉你找不到HHsampnodum.csv

  因为在你的电脑根本没有路径

  C:\Users\Jason\Desktop\Medicare\......

  你必须自行更改路径

dlmread

  Read ASCII-delimited file of numeric data into matrix

  dlmread只读取ASCII file

  所以 HHsampnodum.csv必须改存 HHsampnodum.txt

  才能读到文件

PartDHH4BS.m 有用到function randint()

  这在 Communications System Toolbox

  我没安装Communications System Toolbox

  所以抱歉无法继续

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

parameters,errors,p,q,k,k2,t这些参数如何设定？

你的思惟需要改变

重点在fmincon()及likelihood.m

likelihood.m 是fmincon()的objective function

而fmincon() 的主要目的是求解garch model的parameters

既然parameters都已求出

事后你又把这些parameters,errors,p,q,k,k2,..

带进likelihood function,意义何在

%%%%%%%%%%%%%%%%%%%%%%%%%%%

a=x.boot$std.error^2这一行有什么错误？

x=c(5,4,9,6,21,17,11,20,7,10,21,15,13,16,8)

x.boot=boot(data = x, statistic = function(x, i) var(x), R = 999)

x.boot

names(x.boot)

# [1] "t0"        "t"         "R"         "data"      "seed"      "statistic"

# [7] "sim"       "call"      "stype"     "strata"    "weights"

没std.error,如何取出?

请注意std.error是每次不同的

Bootstrap Statistics :

    original    bias    std. error

t1* 34.31429 -2.465542    7.545234

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

Bootstrap Statistics :

    original    bias    std. error

t1* 34.31429 -2.098413    7.791326

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

Bootstrap Statistics :

    original    bias    std. error

t1* 34.31429 -2.266037    7.826279

std.error=sqrt(var(x.boot$t))

         [,1]

[1,] 7.826279

epoh 发表于 2011-8-7 16:55
Error using ==> dlmread  The file 'C:\Users\Jason\Desktop\Medicare\HHsampnodum.csv'   could not be o ...

epoh老师，您好！
       我想问一下：fmincon()这个函数是matlab自带的还是jplv7中的函数？
如何调用？

em.txt (10.03 KB)

2011-8-8 21:24:31 上传

原程序.txt (5.99 KB)

2011-8-8 21:27:30 上传

epoh老师，您好！
         您看一下下面的程序运行错在什么地方？我尝试修改，还是不能解决问题。
非常感谢！
          附件中有程序。


>     #######################################################
>     ###       R code for Lec 10 Examples                 ###
>     #######################################################
>     
> ### Example 1 (pepper moth)
>
>
>   moth<-function(p,n.obs){
+     n<-sum(n.obs)
+     nc<-n.obs[1]
+     ni<-n.obs[2]
+     nt<-n.obs[3]
+     ntt<-nt
+   
+     cat(p,"\n")
+     pct<-pit<-ptt<-rep(0,20)
+     pct[1]<-p[1]
+     pit[1]<-p[2]
+     ptt[1]<-1-p[1]-p[2]
+    for(i in 2:20){
+     pc.old<-pct[i-1]
+     pi.old<-pit[i-1]
+     pt.old<-ptt[i-1]
+         
+     den<-pc.old^2+2*pc.old*pi.old+2*pc.old*pt.old
+     ncc<-nc*pc.old^2/den
+     nci<-2*nc*pc.old*pi.old/den
+     nct<-2*nc*pc.old*pt.old/den
+     nii<-ni*pi.old^2/(pi.old^2+2*pi.old*pt.old)
+     nit<-2*ni*pi.old*pt.old/(pi.old^2+2*pi.old*pt.old)
+      
+     pct<-(2*ncc+nci+nct)/(2*n)
+     pit<-(2*nii+nit+nci)/(2*n)
+     ptt<-(2*ntt+nct+nit)/(2*n)
+     }
+     return(list(pct=pct,pit=pit,ptt=ptt))
+     }
>  
>    n.obs<-c(85,196,341) # observed data,n_c,n_I,n_T
>    p<-c(1/3,1/3)
>    a<-moth(p,n.obs)
0.3333333 0.3333333
>    pct<-a$pct
>    pit<-a$pit
>    ptt<-a$ptt
>    #convergence diagnostics
>    # statistic R
>    rcc=sqrt( (diff(pct)^2+diff(pit)^2)/(pct[-20]^2+pit[-20]^2) )
>    rcc=c(0,rcc)  #adjusts the length to make the table below
>      
>    d1=(pct[-1]-pct[20])/(pct[-20]-pct[20])
>    d1=c(d1,0)
>    d2=(pit[-1]-pit[20])/(pit[-20]-pit[20])
>    d2=c(d2,0)
>
>    #Table output
>    print(cbind(pct,pit,rcc,d1,d2)[1:9,],digits=5)
           pct     pit        rcc       d1      d2
[1,] 0.333333 0.33333 0.0000e+00 0.042502 0.33659
[2,] 0.081994 0.23741 5.7069e-01 0.036933 0.18765
[3,] 0.071249 0.19787 1.6312e-01 0.036727 0.17780
[4,] 0.070852 0.19036 3.5756e-02 0.036719 0.17624
[5,] 0.070837 0.18902 6.5860e-03 0.036719 0.17595
[6,] 0.070837 0.18879 1.1683e-03 0.036719 0.17589
[7,] 0.070837 0.18875 2.0580e-04 0.036719 0.17588
[8,] 0.070837 0.18874 3.6205e-05 0.036720 0.17587
[9,] 0.070837 0.18874 6.3679e-06 0.036736 0.17587
>
> ### Example 2 Bayes EM: mode of posterior
> #generate observed data
> theta=0.3
> n<-50
> x<-drop(rmultinom(1,n,c((2+theta)/4,(1-theta)/2,theta/4)))
> # prior paramters
> a<-.01
> b<-.01
>
>
> th<-rep(0,20)
> th[1]<-0.2
> for(t in 2:20){
+   num<-x[1]*th[t-1]/(2+th[t-1])+x[3]+a-1
+   den<-x[1]*th[t-1]/(2+th[t-1])+x[2]+x[3]+a+b-2
+   th[t]<-num/den
+   }
> rcc<-sqrt( (diff(th)^2+diff(th)^2)/(th[-20]^2+th[-20]^2) )  
> print(cbind(th,c(0,rcc)),digits=5)
           th           
[1,] 0.20000 0.0000e+00
[2,] 0.19297 3.5148e-02
[3,] 0.18928 1.9131e-02
[4,] 0.18732 1.0366e-02
[5,] 0.18627 5.6016e-03
[6,] 0.18570 3.0227e-03
[7,] 0.18540 1.6298e-03
[8,] 0.18524 8.7843e-04
[9,] 0.18515 4.7333e-04
[10,] 0.18510 2.5502e-04
[11,] 0.18508 1.3739e-04
[12,] 0.18506 7.4013e-05
[13,] 0.18506 3.9871e-05
[14,] 0.18505 2.1478e-05
[15,] 0.18505 1.1570e-05
[16,] 0.18505 6.2329e-06
[17,] 0.18505 3.3576e-06
[18,] 0.18505 1.8087e-06
[19,] 0.18505 9.7435e-07
[20,] 0.18505 5.2488e-07
> ### Example 3 SEM algorithm
>
> #First, run EM using the code in the basic EM algorithm given in Example 1.
> #Then we run SEM as follows.  
>
> #Starting values for SEM
> #Choosing the parameter values at the 20th iteration of EM.  
> pcthat=pct[20]   
> pithat=pit[20]
> ptthat=ptt[20]
>
> #Initialize the parameter estimate vectors, one vector for each phenotype
> sempct=rep(0,20)  
> sempit=rep(0,20)  
> semptt=rep(0,20)
> sempct[1]=0.07      #Initial SEM estimates of the parameters
> sempit[1]=0.19
> semptt[1]=0.74
> rij=array(0,c(3,3,20))
>
> #This for loop implements the SEM algorithm for the peppered moth
> #example.  
>
> #Below t is the number of SEM iterations
> for (t in 2:20) {
+   #Take standard EM step
+   denom1=(sempct[t-1]^2+2*sempct[t-1]*sempit[t-1]+2*sempct[t-1]*semptt[t-1])
+   ncct=nc*sempct[t-1]^2/denom1
+   ncit=2*nc*sempct[t-1]*sempit[t-1]/denom1
+   nctt=nc-ncct-ncit
+   niit=ni*sempit[t-1]^2/(sempit[t-1]^2+2*sempit[t-1]*semptt[t-1])
+   nitt=ni-niit
+   nttt=nt
+   sempct[t]=(2*ncct+ncit+nctt)/(2*n)
+   sempit[t]=(2*niit+ncit+nitt)/(2*n)
+   semptt[t]=(2*nttt+nctt+nitt)/(2*n)
+   #SEM
+   #loop over each parameter
+   for (j in 1:3) {
+     #start at estimates from the the 20th iteration of EM
+     sempj=c(pcthat,pithat,ptthat)  
+
+     #replace the jth element of sempj with the most recent EM estimate
+     sempj[j]=c(sempct[t],sempit[t],semptt[t])[j]
+
+     #Take one EM step for sempj   
+     denom1=(sempj[1]^2+2*sempj[1]*sempj[2]+2*sempj[1]*sempj[3])
+     ncct=nc*sempj[1]^2/denom1
+     ncit=2*nc*sempj[1]*sempj[2]/denom1
+     nctt=nc-ncct-ncit
+     niit=ni*sempj[2]^2/(sempj[2]^2+2*sempj[2]*sempj[3])
+     nitt=ni-niit
+     nttt=nt
+     nextstep=c((2*ncct+ncit+nctt)/(2*n),(2*niit+ncit+nitt)/(2*n),
+        (2*nttt+nctt+nitt)/(2*n))
+
+     # Calculate rij.  This is (4.44) on page 101
+     rij[,j,t]=(nextstep-c(pcthat,pithat,ptthat))/
+                   (sempj[j]-c(pcthat,pithat,ptthat)[j])
+       }
+   }
错误: 找不到对象'nc'
>   for(t in 1:20){
+    cat(t,sempct[t],sempit[t],semptt[t])
+    cat("\n")
+    print(rij[,,t])
+    cat("\n")
+    }
1 0.07 0.19 0.74
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0    0

2 0 0 0
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0    0

3 0 0 0
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0    0

4 0 0 0
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0    0

5 0 0 0
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0    0

6 0 0 0
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0    0

7 0 0 0
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0    0

8 0 0 0
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0    0

9 0 0 0
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0    0

10 0 0 0
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0    0

11 0 0 0
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0    0

12 0 0 0
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0    0

13 0 0 0
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0    0

14 0 0 0
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0    0

15 0 0 0
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0    0

16 0 0 0
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0    0

17 0 0 0
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0    0

18 0 0 0
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0    0

19 0 0 0
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0    0

20 0 0 0
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0    0

> #Note that the algorithm becomes unstable on 8th iteration
> #Below is the output for iteration 7
> #EM after 20 iterations
> #> cbind(pct,pit,ptt)[20,]
> #       pct        pit        ptt
> #0.07083691 0.18873652 0.74042657
>
> #SEM after 7 iterations (the way the indices work, this is index 6
> #> cbind(sempct,sempit,semptt)[6,]
> #    sempct     sempit     semptt
> #0.07083691 0.18873670 0.74042639
>
> rij[,,6]
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0    0
> #             [,1]         [,2]        [,3]
> #[1,]  0.034117920 -0.002601059 -0.00260106
> #[2,] -0.006910223  0.140676982 -0.03519589
> #[3,] -0.027207449 -0.138075923  0.03779695
>
> #Now need iyhat (inverse of the information matrix for Y)  used in (4.43)
>
> #Equations (4.5)-(4.9) on page 93
> denom1=(pcthat^2+2*pcthat*pithat+2*pcthat*ptthat)
> ncct=nc*pcthat^2/denom1
错误: 找不到对象'nc'
> ncit=2*nc*pcthat*pithat/denom1
错误: 找不到对象'nc'
> nctt=nc-ncct-ncit
错误: 找不到对象'nc'
> niit=ni*pithat^2/(pithat^2+2*pithat*ptthat)
错误: 找不到对象'ni'
> nitt=ni-niit
错误: 找不到对象'ni'
> nttt=nt
错误: 找不到对象'nt'
>

> #Required derivatives for iyhat (used in eqn (4.41) on page 101)
> d20q=-(2*ncct+ncit+nctt)/pcthat^2 -(2*nttt+nitt+nctt)/(ptthat^2)
错误: 找不到对象'ncct'
> d02q=-(2*niit+ncit+nitt)/pithat^2 -(2*nttt+nitt+nctt)/(ptthat^2)
错误: 找不到对象'niit'
> d12q=-(2*nttt+nitt+nctt)/(ptthat^2)
错误: 找不到对象'nttt'
>  
> iyhat=-cbind(c(d20q,d12q),c(d12q,d02q))
错误于cbind(c(d20q, d12q), c(d12q, d02q)) : 找不到对象'd20q'
>
> solve(iyhat)
错误于solve(iyhat) : 找不到对象'iyhat'
> #[1,]  5.290920e-05 -1.074720e-05
> #[2,] -1.074720e-05  1.230828e-04
>
> #Since the percentages of the 3 moth phenotypes add up to 1, we only
> #presented results for carbonaria and insularia phenotypes.
>
> psiprime=rij[,,6] #(7th iteration)
> psiprime22=psiprime[-3,-3]
>
> #variance matarix of MLE
> varhat=solve(iyhat)%*%solve(diag(2)-t(psiprime22))
错误于solve(iyhat) : 找不到对象'iyhat'
> varhat=(varhat+t(varhat))/2
错误: 找不到对象'varhat'
> varhat
错误: 找不到对象'varhat'
> #              [,1]          [,2]
> #[1,]  5.481298e-05 -1.223007e-05
> #[2,] -1.223007e-05  1.433249e-04
>
> #sd(pc), sd(pi)
> sqrt(diag(varhat))              
错误于diag(varhat) : 找不到对象'varhat'
> #[1] 0.007403579 0.011971838
>
>
> #cor(pc,pi)
> varhat[1,2]/prod(sqrt(diag(varhat)))  
错误: 找不到对象'varhat'
> #[1] -0.1379833
>
> #var(pt)
> varhat[1,1]+varhat[2,2]+2*varhat[1,2]  
错误: 找不到对象'varhat'
> #[1] 0.0001736777
>
> #sd(pt)
> sqrt(sum(varhat))
错误: 找不到对象'varhat'
> #[1] 0.01317868
>
> #cor(pc,pt)
> (-varhat[1,1]-varhat[1,2])/(sqrt(varhat[1,1])*sqrt(sum(varhat)))
错误: 找不到对象'varhat'
> #[1] -0.4364370
>
> #cor(pi,pt)
> (-varhat[2,2]-varhat[1,2])/(sqrt(varhat[2,2])*sqrt(sum(varhat)))
错误: 找不到对象'varhat'
> #[1] -0.8309075
>

EGARCH.txt (88.66 KB)

2011-8-8 21:38:37 上传

epoh老师，您好！
        您看一下附件中的EGARCH程序是不是用jplv7中的ucsd-garch工具箱可以运行？
我运行总是出错，您看如何运行？
非常感谢！

怎么也没找到Gauss的code,但是找到了一个Matlab code.是一篇德国人写的论文，现介绍了模型，然后还应用了德国的数据，最后提供了Matlab code。大家看看吧~真遗憾，Gauss 竟然没有现成的程序。
Markov-switching EGARCH (Matlab).pdf

本文来自: 人大经济论坛 Gauss专版 版，详细出处参考： https://bbs.pinggu.org/forum.php? ... 1&from^^uid=11232