copula的极大似然估计

请教高手啊，copula 的极大似然估计怎么用matlab实现？

1. u=data(:,1);
   
2. v=data(:,2);
   
3. T=length(u);
   
4. options = optimset('Algorithm','interior-point','Display','iter','TolCon',10^-12,'TolFun',10^-4,'TolX',10^-6); %优化条件
   
5. 
   
6. 
   
7. % 1. Normal Copula
   
8. kappa1 = corrcoef12(norminv(u),norminv(v));
   
9. LL1 = NormalCopula_CL(kappa1,[u,v]);          
   
10. 
    
11. % 2. Clayton's copula
    
12. lower = 0.0001;
    
13. theta0 = 1;
    
14. [ kappa2 LL2] = fmincon('claytonCL',theta0,[],[],[],[],lower,[],[],options,[u,v]); %参数约束非线性优化方法估计clayton copula函数的参数
    
15. 
    
16. % 3. Rotated Clayton copula (with tail dep in upper tail instead of lower)
    
17. lower = 0.0001;
    
18. theta0 = 1;
    
19. [ kappa3 LL3] = fmincon('claytonCL',theta0,[],[],[],[],lower,[],[],options,1-[u,v]);
    
20. 
    
21. % 4. Plackett copula
    
22. lower = 0.0001;
    
23. theta0 = 1;
    
24. [ kappa4 LL4] = fmincon('plackettCL',theta0,[],[],[],[],lower,[],[],options,[u,v]);
    
25. % LL5 = -3.2721
    
26. 
    
27. % 5. Frank copula
    
28. theta0 = 1;
    
29. [ kappa5 LL5] = fmincon('frankCL',theta0,[],[],[],[],lower,[],[],options,[u,v]);
    
30. 
    
31. % 6. Gumbel copula
    
32. lower = 1.1;
    
33. theta0 = 2;
    
34. [ kappa6 LL6] = fmincon('gumbelCL',theta0,[],[],[],[],lower,[],[],options,[u,v]);
    
35. 
    
36. % 7. Rotated Gumbel copula
    
37. lower = 1.1;
    
38. theta0 = 2;
    
39. [ kappa7 LL7] = fmincon('gumbelCL',theta0,[],[],[],[],lower,[],[],options,1-[u,v]);
    
40. 
    
41. % 8. Student's t copula
    
42. lower = [-0.999 , 2.1 ];
    
43. upper = [ 0.999 , 100 ];
    
44. theta0 = [kappa1;10];
    
45. [ kappa8 LL8] = fmincon('tcopulaCL',theta0,[],[],[],[],lower,upper,[],options,[u,v]);
    
46. 
    
47. % 9. Symmetrised Joe-Clayton copula
    
48. lower = [0 , 0 ];
    
49. upper = [ 1 , 1];
    
50. theta0 = [0.25;0.25];
    
51. [ kappa9 LL9] = fmincon('sym_jc_CL',theta0,[],[],[],[],lower,upper,[],options,[u,v]);
    
52. 
    
53. LL = [LL1;LL2;LL3;LL4;LL5;LL6;LL7;LL8;LL9];

复制代码

u,v为服从0,1均匀分布序列，'tcopulaCL'等等为编写好copula似然函数

怎么没有回复的啊 。。。。

[程序问答] 请高手看看，matlab 极大似然估计程序哪出错了 [复制链接]  
  

tiaotiaotang tiaotiaotang 当前离线
当日剩余下载次数23 次阅读权限20论坛币7 个经验330 点威望0 级帖子47精华0学术水平1 点热心指数1 点信用等级1 点在线时间143 小时注册时间2007-12-14最后登录2011-12-12
加关注
0
关注 0
粉丝 高中生



高中生, 积分 47, 距离下一级还需 3 积分论坛币7 个经验330 点威望0 级帖子47精华0学术水平1 点热心指数1 点信用等级1 点在线时间143 小时注册时间2007-12-14最后登录2011-12-12串个门 加好友 打招呼 发消息  电梯直达   楼主
发表于 昨天 10:56 |只看该作者 |倒序浏览 |推荐到门户首页
我做的过程是：
将y.xls导入matlab
将skewtdis_LL.m, my_mle.m加入MATLAB.然后在命令窗口输入
[para,standard_deviation]=my_mle('skewtdis_LL',[5;0.02],y)
结果显示：
??? Error using ==> skewtdis_LL
Too many input arguments.
Error in ==> fminsearch at 195
fv(:,1) = funfcn(x,varargin{:});
Error in ==> my_mle at 3
[para,fv]=fminsearch(fun,para0,[],2,varargin{:});


程序如下：
function LL = skewtdis_LL(theta,x)
% PURPOSE: returns the log-likelihood at x of Hansen's (1994) 'skewed t' distribution
%---------------------------------------------------
% USAGE: LL = skewtdis_LL(x,nu,lambda)
% where: x  = a vector of data
%        theta = [nu;lambda]
%     nu = degrees of freedom parameter
%       lambda = skewness parameter
%---------------------------------------------------
% RETURNS:
%        a matrix of log-likelihood at each element of x
% --------------------------------------------------
% SEE ALSO: tdis_cdf, tdis_rnd, tdis_inv, tdis_prb, skewtdis_pdf
%---------------------------------------------------
%
%  Andrew Patton
%
%  25 June, 2001
% This code was used in:
%
% Patton, Andrew J., 2002, "On the Out-of-Sample
% Importance of Skewness and Asymmetric Dependence for
% Asset Allocation", working paper, Department of Economics,
% University of California, San Diego.
nu = theta(1);
lambda = theta(2);
[T,k] = size(x);
nu = nu(1)*ones(T,1);     % can make this time-varying, but needs to be >2
lambda = lambda(1)*ones(T,1); % can make this time-varying, but needs to be in (-1,1)
   
c = gamma((nu+1)/2)./(sqrt(pi*(nu-2)).*gamma(nu/2));
a = 4*lambda.*c.*((nu-2)./(nu-1));
b = sqrt(1 + 3*lambda.^2 - a.^2);
logc = gammaln((nu+1)/2) - gammaln(nu/2) - 0.5*log(pi*(nu-2));
logb = 0.5*log(1 + 3*lambda.^2 - a.^2);
find1 = (x<(-a./b));
find2 = (x>=(-a./b));
LL1   = logb + logc - (nu+1)/2.*log(1+1./(nu-2).*((b.*x+a)./(1-lambda)).^2);
LL2   = logb + logc - (nu+1)/2.*log(1+1./(nu-2).*((b.*x+a)./(1+lambda)).^2);
LL    = sum(LL1(find1)) + sum(LL2(find2));
LL    = -LL;
function [para,standard_deviation,fv]=my_mle(fun,para0,varargin)
para0=para0(:);
[para,fv]=fminsearch(fun,para0,[],2,varargin{:});
fv=-fv;
d=numericalfirstderivative(fun,para,1,varargin{:});
standard_deviation=sqrt(diag(pinv(d'*d)));
function f=numericalfirstderivative(fun,parameter,varargin)
n=length(parameter);
for i=1:n
a=zeros(n,1);
a(i)=min(parameter(i)*1e-6,1e-5);
y1(:,i)=feval_r(fun,parameter+a,varargin{:});
y2(:,i)=feval_r(fun,parameter-a,varargin{:});
f(:,i)=(y1(:,i)-y2(:,i))/2/a(i);
end



淘宝网购物自动返论坛币  
skewtdis_LL.doc
昨天 10:56 上传

21.5 KB


回帖推荐epoh 发表于2楼  查看完整内容

% skewtdis_rnd.m need Ucsd_garch toolbox x=skewtdis_rnd(8.1,-0.09,10000); init=[8.5;-0.2] lb=[0,-1] ub= =fmincon(@(theta)skewtdis_LL(theta,x),init,[],[],[],[],lb,ub) para = 8.3054 -0.0837
分享到: QQ空间 腾讯微博 腾讯朋友 收藏0 评分  
回复 使用道具 举报
提升卡

变色卡
  
  

epoh epoh 当前离线
当日剩余下载次数29 次阅读权限30论坛币4908 个经验3574 点威望1 级帖子739精华0学术水平1478 点热心指数1483 点信用等级1432 点在线时间2 小时注册时间2006-11-26最后登录2011-12-12
加关注
0
关注 24
粉丝 副教授



副教授, 积分 739, 距离下一级还需 151 积分论坛币4908 个经验3574 点威望1 级帖子739精华0学术水平1478 点热心指数1483 点信用等级1432 点在线时间2 小时注册时间2006-11-26最后登录2011-12-12串个门 加好友 打招呼 发消息  沙发
发表于 昨天 16:05 |只看该作者
% skewtdis_rnd.m need Ucsd_garch toolbox
x=skewtdis_rnd(8.1,-0.09,10000);
init=[8.5;-0.2]
lb=[0,-1]
ub=[Inf,0]

[para,llf,exitflag]=fmincon(@(theta)skewtdis_LL(theta,x),init,[],[],[],[],lb,ub)
para =

    8.3054
   -0.0837

已有 3 人评分 学术水平 热心指数 信用等级 收起 理由  
zhangtao  + 5 + 5 + 5 精彩帖子
tiaotiaotang  + 1 + 1  热心，水平高
liuxin9023  + 1   分析的有道理
总评分: 学术水平 + 7  热心指数 + 6  信用等级 + 5   查看全部评分


回复  评分 举报

  

tiaotiaotang tiaotiaotang 当前离线
当日剩余下载次数23 次阅读权限20论坛币7 个经验330 点威望0 级帖子47精华0学术水平1 点热心指数1 点信用等级1 点在线时间143 小时注册时间2007-12-14最后登录2011-12-12
加关注
0
关注 0
粉丝 高中生



高中生, 积分 47, 距离下一级还需 3 积分论坛币7 个经验330 点威望0 级帖子47精华0学术水平1 点热心指数1 点信用等级1 点在线时间143 小时注册时间2007-12-14最后登录2011-12-12串个门 加好友 打招呼 发消息  藤椅
发表于 昨天 22:37 |只看该作者
epoh 发表于 2011-12-11 16:05
% skewtdis_rnd.m need Ucsd_garch toolbox
x=skewtdis_rnd(8.1,-0.09,10000);
init=[8.5;-0.2]太感谢了，解决我的大问题了


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

谢了。。

同求

yucongy 发表于 2011-11-24 21:43
u,v为服从0,1均匀分布序列，'tcopulaCL'等等为编写好copula似然函数

你好 ，还在弄'tcopulaCL'等等为编写好copula似然函数相关文章么，我也在弄。我的程序运行到

41.% 8. Student's t copula

42.lower = [-0.999 , 2.1 ];

43.upper = [ 0.999 , 100 ];

44.theta0 = [kappa1;10];

45.[ kappa8 LL8] = fmincon('tcopulaCL',theta0,[],[],[],[],lower,[],[],options,1-[u,v]);


之后就提示出错了，不知道错在哪里？谢谢额，希望给与指点！

sisters 发表于 2013-1-16 16:21
你好 ，还在弄'tcopulaCL'等等为编写好copula似然函数相关文章么，我也在弄。我的程序运行到

41.% 8.  ...

写上错误的提示说明，这样才方便解答

yucongy 发表于 2013-4-21 13:43
写上错误的提示说明，这样才方便解答

Your initial point x0 is not between bounds lb and ub; FMINCON
shifted x0 to satisfy the bounds.
Undefined function 'tdis_inv' for input arguments of type 'double'.

Error in tcopulaCL (line 33)
   x = tdis_inv(data(:,1),nu);

Error in fmincon (line 631)
      initVals.f = feval(funfcn{3},X,varargin{:});

Caused by:
    Failure in initial user-supplied objective function evaluation. FMINCON cannot contin
我也遇到同样的问题，错误提示是上面所写的，请问怎么解决？谢谢~

非常感谢分享