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

88# epoh


# MNZModel.ssc                        Morley, Nelson and Zivot's UC model
#
# Created: January 8, 2005

##
## MNZ model: Ref. Morley, Nelson and Zivot (2002), "Why are
## Beveridge-Nelson and Unobserved Components Decompositions of GDP
## so Different?", forthcoming in Review of Economics and Statistics
##

## Note: parameters are taken from Table 4 of
## Morley, Nelson and Zivot (2002)
##
delta = 0.8156
phi1 = 1.3419
phi2 = -0.7060
sigma.v = 1.2368
sigma.w = 0.7485
sigma.vw = -0.8389
bigV = matrix(c(sigma.v^2,sigma.vw,sigma.vw,sigma.w^2),2,2)
Omega = matrix(0,4,4)
Omega[1:2,1:2] = bigV
a1 = matrix(0,3,1)
# solve for initial variance of stationary part
bigF = matrix(c(phi1,1,phi2,0),2,2)
vecV = c(sigma.w^2,0,0,0)
vecP = solve(diag(4)-kronecker(bigF,bigF))%*%vecV
P.ar2 = matrix(vecP,2,2)
Sigma = matrix(0,4,3)
Sigma[1,1] = -1
Sigma[2:3,2:3] = P.ar2
ssf.mnz= list(mDelta=c(delta,0,0,0),
mPhi=rbind(c(1,0,0),c(0,phi1,phi2),c(0,1,0),c(1,1,0)),
mOmega=Omega,
mSigma = Sigma)
ssf.mnz

## simulate from model
##
set.seed(569)
mnz.sim = SsfSim(ssf.mnz,n=250)
par(mfrow=c(2,1))
tsplot(mnz.sim[,c(1,4)],main="Simulated response and trend")
tsplot(mnz.sim[,2],main="Simulated cycle")

## estimate model using US postwar quarterly real GDP
##
MNZ.mod = function(parm) {
        delta = parm[1]
        phi1 = parm[2]
        phi2 = parm[3]
        sigma.v = exp(parm[4])
        sigma.w = exp(parm[5])
        rho.vw = parm[6]
        sigma.vw = sigma.v*sigma.w*rho.vw
        bigV = matrix(c(sigma.v^2,sigma.vw,sigma.vw,sigma.w^2),2,2)
        Omega = matrix(0,4,4)
        Omega[1:2,1:2] = bigV
        a1 = matrix(0,3,1)
# solve for initial variance of stationary part
        bigF = matrix(c(phi1,1,phi2,0),2,2)
        vecV = c(sigma.w^2,0,0,0)
        vecP = solve(diag(4)-kronecker(bigF,bigF))%*%vecV
        P.ar2 = matrix(vecP,2,2)
        Sigma = matrix(0,4,3)
        Sigma[1,1] = -1
        Sigma[2:3,2:3] = P.ar2
        ssf.mod= list(mDelta=c(delta,0,0,0),
        mPhi=rbind(c(1,0,0),c(0,phi1,phi2),c(0,1,0),c(1,1,0)),
        mOmega=Omega,
        mSigma = Sigma)
        CheckSsf(ssf.mod)
}

## set starting values
##
MNZ.start=c(0.81,1.34,-0.70,0.21,-0.30,-0.9)
names(MNZ.start) = c("delta","phi.1","phi.2",
"ln.sigma.v","ln.sigma.w","rho")

## estimate model by MLE
##
low.vals = c(0,0,-2,-Inf,-Inf,-0.999)
up.vals = c(2,2,2,Inf,Inf,0.999)
MNZ.mle = SsfFit(MNZ.start,lny.ts,"MNZ.mod",
lower=low.vals,upper=up.vals)

## estimated parameters and asymptotic standard errors
##
MNZ.mle$parameters
sqrt(diag(MNZ.mle$vcov))
# estimated standard deviations
exp(MNZ.mle$parameters[4:5])
# note: use delta method to get standard errors for standard deviations

# compute filtered trend and cycle
filteredEst.MNZ = SsfMomentEst(lny.ts,
MNZ.mod(MNZ.mle$parameters),task="STFIL")
trend.filter = timeSeries(data=filteredEst.MNZ$state.moment[,1],
positions=filteredEst.MNZ$positions)
cycle.filter = timeSeries(data=filteredEst.MNZ$state.moment[,2],
positions=filteredEst.MNZ$positions)

par(mfrow=c(2,1))
plot(trend.filter,lny.ts,
main="Log Real GDP and Filtered Trend from MNZ Model",reference.grid=F)
plot(cycle.filter,main="Filtered Cycle from MNZ Model",reference.grid=F)

# compute smoothed trend and cycle
smoothedEst.MNZ = SsfMomentEst(lny.ts,
MNZ.mod(MNZ.mle$parameters),task="STSMO")
trend.smooth = timeSeries(data=smoothedEst.MNZ$state.moment[,1],
positions=smoothedEst.MNZ$positions)
cycle.smooth = timeSeries(data=smoothedEst.MNZ$state.moment[,2],
positions=smoothedEst.MNZ$positions)

par(mfrow=c(2,1))
plot(trend.smooth,lny.ts,
main="Log Real GDP and Smoothed Trend from MNZ Model")
plot(cycle.smooth,main="Smoothed Cycle from MNZ Model")

plot(cycle.smooth,reference.grid=F)

# plot smoothed cycle with 95% error bands
cycle.sd= sqrt(timeSeries(data=smoothedEst.MNZ$state.variance[,2],
positions=smoothedEst.MNZ$positions))
upper = cycle.smooth + 2*cycle.sd
lower = cycle.smooth - 2*cycle.sd
par(mfrow=c(1,1))
plot(cycle.smooth,upper,lower,
plot.args=list(lty=c(1,3,3)))

# plot filtered and smoothed cycle together
plot(cycle.filter,cycle.smooth,
plot.args=list(lty=c(1,2)))
legend(0.75,-3,legend=c("filtered","smoothed"),lty=c(1,2))

# plot filtered and smoothed trend together
plot(diff(trend.filter),diff(trend.smooth))

今天系统附件传不来程序，所以粘出来，估计这个程序和数据epoh老师您有
运行以上程序，为什么会提示以下错误：
> mnz.sim = SsfSim(ssf.mnz, n = 250)
Problem: Couldn't find a function definition for "SsfSim"
Use traceback() to see the call stack

##
low.vals = c(0, 0, -2,  - Inf,  - Inf, -0.999)
> up.vals = c(2, 2, 2, Inf, Inf, 0.999)
> MNZ.mle = SsfFit(MNZ.start, lny.ts, "MNZ.mod", lower = low.vals,
        upper = up.vals)
Problem: Object "lny.ts" not found
Use traceback() to see the call stack
装入Finmetrics模块后，又提示以上错误，我没有办法了，请epoh老师帮忙！
非常感谢！

traceback() 这个函数有什么用？如何用？
> traceback()
14: eval(action, sys.parent())
13: doErrorAction("Problem: Object \"lny.ts\" no
t found", 1000)
12: as.character(class(x))
11: match(el, set, nomatch = NA)
10: is.na(match(el, set, nomatch = NA))
9: !is.na(match(el, set, nomatch = NA))
8: is.element(as.character(class(x)), c("bdChara
cter", "bdExclude",
7: stop.on.bdObject(data)
6: SsfFit(MNZ.start, lny.ts, "MNZ.mod", lower =
low.vals, upper =
5: eval(i, local)
4: source(auto.print = auto.print, exprs = subst
itute(exprs.literal))
3: script.run(exprs.literal = {
2: eval(expression(script.run(exprs.literal = {
1:
Message: Problem: Object "lny.ts" not found
>
以上结果如何看和分析？

缺少数据lngdpq.txt
lngdpq=read.table("lngdpq.txt")
lngdpq.df=lngdpq[,1]

lny = as.matrix(lngdpq.df)*100

td = timeSeq(from="1/1/1947",to="4/1/1998",by="quarters")
lny.ts = timeSeries(data=lny,positions=td)
source("MNZUCModel.ssc")

lngdpq.txt

   lngdpq.txt (2.79 KB)

2011-7-25 16:14:54 上传

###########
当你要执行checkSV.ssc
首先要先打开checkSV.ssc
看要用到哪些数据及函数
State Space Modeling
把这些摆在一起时

StateSpaceSurvey.ssc; sv_mcl_est.ssc; ssfnong.ssc; checkSV.ssc
他们之间大都会相互调用函数.

94# epoh

经您指导后，运行非常好！非常感谢！
lny = as.matrix(lngdpq.df)*100
我想问一下：为什么要在上面一行乘以100，然后再做成一个矩阵？

td = timeSeq(from="1/1/1947",to="4/1/1998",by="quarters")
lny.ts = timeSeries(data=lny,positions=td)
为什么在上面两行要对数据指定位置和做成时间序列，难道原来的数据不是时间序列吗？


from="1/1/1947",to="4/1/1998",by="quarters"
以上数据的开始时间和最后时间，季节数据您是怎么知道的？我在原始数据和程序中怎么没有看出来？



缺少数据lngdpq.txt
lngdpq=read.table("lngdpq.txt")
lngdpq.df=lngdpq[,1]

lny = as.matrix(lngdpq.df)*100

td = timeSeq(from="1/1/1947",to="4/1/1998",by="quarters")
lny.ts = timeSeries(data=lny,positions=td)
source("MNZUCModel.ssc")

lngdpq.txt

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

StateSpaceSurvey.ssc;MNZUCModel.ssc; BNdecomposition.ssc;
都有用到lny.ts

lny = as.matrix(lngdpq.df)*100        
td = timeSeq(from="1/1/1947",to="4/1/1998",by="quarters")
lny.ts = timeSeries(data=lny,positions=td)
以上三行取自BNdecomposition.ssc

log real GDP data:lngdpq.txt
取自
  http://faculty.washington.edu/ezivot/MFTS2ndEditionAdditions.htm


ps:

data=read.table()
data 是一个"data.frame"
由于会用到tsplot画图
所以用timeSeries()转换为"timeSeries" object


通常optimization过程,为了顺利收敛.
都会100*data
就像估计garch model
都采100* logreturn

epoh老师，您好！
    我执行：markovSwitchingExamples.SSC时提示以下错误：Problem in parse(text = "script.run(exprs.literal={\n## markovSwitc..: Syntax error: Unbalanced parentheses, expected ")" before unbalanced ("}")  on input line 148
Use traceback() to see the call stack
想尽了一切办法无法解决，您看如何解决呢？
非常感谢！

65# epoh
请把full_bekk_mvgarch.m
这三行:
scalaropt=optimset('fminunc');
scalaropt=optimset(scalaropt,'TolFun',1e-1,'Display','iter','Diagnostics','on','DiffMaxChange',1e-2);
startingparameters=scalar_bekk_mvgarch(data,p,q,scalaropt);

改为这四行:

%scalaropt=optimset('fminunc');
%scalaropt=optimset(scalaropt,'TolFun',1e-1,'Display','iter','Diagnostics','on','DiffMaxChange',1e-2);
%startingparameters=scalar_bekk_mvgarch(data,p,q,scalaropt);

尊敬的epoh老师，您好！
    为什么以上三行注释了以后，程序就可以正常运行了呢？

startingparameters=[0.6836; 0.2031; 0.3717; 0.2803; 0.9101]
   
为什么要对startingparameters设置初值，而且必须要设五个[0.6836; 0.2031; 0.3717; 0.2803; 0.9101]，而不是四个初值：[0.6836; 0.2031; 0.3717; 0.2803]？
非常感谢！

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

96# epoh

epoh老师，还有附件中的这个VaR程序，运行后提示以下错误，我实在没办法，还希望
老师您帮忙看看！
非常感谢！
Problem in parse(text = "script.run(exprs.liter..: Syntax error: Unbalanced parentheses, expected "}" before unbalanced (")")  on input line 126
Use traceback() to see the call stack

###"markovSwitchingExamples.ssc"
1.漏掉一个 ")",我已补上.

2.少了数据ndx.rvol我也已补上.
  ndx.ret2 = getReturns(ndx.dat[,"Close"])^2
  ndx.rvol = sqrt(aggregate(ndx.ret2, FUN=sum, by="weeks", week.align=1))

3.少了function shadePlot()
  遍寻不着,就在程序最后面.
  我把他封了.
   markovSwitchingExam.rar (3.04 KB) 本附件包括：

• markovSwitchingExamples.SSC

2011-7-26 15:10:20 上传

###########"valueAtRiskDemo.ssc"

错在底下norm.ETL多了一个")"


norm.ETL = function(x, VaR) {
## required arguments:
## x   timeSeries of returns (simple or continuous)
## VaR   VaR estimate based on normal approximation.
## optional arguments:
## value:
##    numeric value giving ETL estimate from estimated normal
##   distribution
mu = colMeans(x)
sd = colStdevs(x)
z = (mu - VaR)/sd)
ETL = mu + sd*(dnorm(z)/pnorm(z))
ETL
}

请将z = (mu - VaR)/sd)

改为z = (mu - VaR)/sd



%%%%%%%%%%matlab full_bekk
正常的话,startingparameters应该由scalar_bekk_mvgarch估计
但因为当时你的function vech()出问题
所以才有这个临时的权宜之计
先供你能跑程序
不知你重安装之后,是否正常了?

full_bekk     要估计的系数有11个
diagonal_bekk 要估计的系数有 9个
scalar_bekk   要估计的系数有 5个

以full_bekk 而言
startingparameters=[0.6836; 0.2031; 0.3717; 0.2803; 0.9101]
只是第一步而已,因为是要估计11个系数
所以startingparameters也是11个
startingparameters=[0.6836; 0.2031; 0.3717; 0.2803; 0;0;0.2803;0.9101;0;0;0.9101]
透过full_bekk_mvgarch_likelihood.m 来估计
parameters=fminunc('full_bekk_mvgarch_likelihood',startingparameters,options,data,p,q,k,k2,t);