半参数以及混合藤copula 的MATLAB程序分享

相似文件 换一批

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

已经下载论坛很多书籍，大概6000论坛币已经剩下只有300左右吧，此贴不收金币。我太累，已经看不了那么多书。

tulipsliu  中国数量经济（群）  中国数量金融（群）  群主

原则上中国数量经济已经封闭群，只实行邀请制，没通过我的邀请无法进群。

中国数量金融群为开放群，这里分享群号码，想一起学习金融的可以加入。

群嘛，不知道怎么介绍。中国数量金融群号码：949614744

其中应该是半参数那个，作者在程序里设置了  Frank copula ,程序运行中途无法继续。我已经修改为 t copula
可以放心运行。贴上半参数里其中一个主程序的代码，可以不复制到你的文件里，因为程序全部在 ZIP 文件里。
这个程序应该是我在 github 网站下载，如果想看更新最新的代码，可以自己去搜索。

1. %This MATLAB script runs the Semi-Parametric Copula Method presented in Manheim &
   
2. %Detwiler (2019) for various analytical test cases. It was originally coded and developed by Derek Manheim,
   
3. %at the University of California, Irvine, Department of Civil & Environmental
   
4. %Engineering using MATLAB r2015b. It is a fully Global, Moment Independent, Monte Carlo (MC) based sampling method
   
5. %that relies on the use of parametric copula models to describe the joint pdf between
   
6. %model inputs (i.e., parameters) and outputs. The user must specify what
   
7. %test case to run (testfunc, options are 1,2,3), the copula model used
   
8. %(options are 'Gaussian','t','Clayton','Frank','Gumbel'), the number of MC sampling points used to construct
   
9. %and sample the Copula model (Nevalpts) as well as the repetition number (Rep), which specifies the
   
10. %location of the Sobol sequence used for QMC sampling. In addition to the files provided in this download folder,
    
11. %in order to run this code users must have the statistics and machine learning toolbox package installed
    
12. %and should be running MATLAB r2015 or later. The results are automatically
    
13. %saved after running this script. *Note - running over 65536 Nevalpts will
    
14. %require more advanced parallel computation (8-16 cores or more) and is not recommended to be run locally.
    
15. %Final Version, updated 5-8-2019
    
16. clear ; close all; clc
    
17. 
    
18. %Run specifications
    
19. Rep = 1; %Location in the sobol sequence (1,2,3,4) to set up different repetitions
    
20. testfunc = 1; %Which analytical test cases?
    
21. CopulaModel = 't'; %Which parametric Copula Model to use?
    
22. Nevalpts = 256; %How many MC evaluation points used to construct and sample the Copula model? should be a factor of 2
    
23. Delay=500+((Rep-1)*Nevalpts); %This sets the location in the sobol sequence for QMC sampling
    
24. 
    
25. %Set up variables required for MLE optimization
    
26. if strcmp(CopulaModel,'Clayton')
    
27.     model = 1;
    
28. elseif strcmp(CopulaModel,'Frank')
    
29.     model = 2;
    
30. elseif strcmp(CopulaModel,'Gumbel')
    
31.     model = 3;
    
32. elseif strcmp(CopulaModel,'Gaussian')
    
33.     model = 4;
    
34. elseif strcmp(CopulaModel,'t')
    
35.     model = 5;
    
36. end
    
37. 
    
38. %%Set up variables to evaluate each test case
    
39. if testfunc == 1 %linear test case
    
40.    Nvar = 6;
    
41.    f = @(Xs) (1.5.*Xs(:,1)+1.6.*(Xs(:,2))+1.7.*(Xs(:,3))+1.8.*(Xs(:,4))+1.9.*(Xs(:,5))+2.*(Xs(:,6)));
    
42.    p = sobolset(Nvar,'skip',Delay);
    
43.    X = p(1:Nevalpts,:);
    
44.    for j = 1:Nvar
    
45.        Xs(:,j) = sqrt(2)*erfinv(2.*X(:,j)-1); %Model inputs
    
46.    end
    
47.   [Yn] = f(Xs);
    
48.    Y = Yn; %Model output
    
49. elseif testfunc == 2 %non-linear test case
    
50.     Nvar = 6;
    
51.     bi = [4.5,2.5,1.5,0.5,0.3,0.1];
    
52.     f = @(x,bi) exp(sum(bi.*x,2))-prod((exp(bi)-1)./bi);
    
53.     p = sobolset(Nvar,'skip',Delay);
    
54.     X = p(1:Nevalpts,:);
    
55.     for j = 1:Nevalpts
    
56.         Yn(j,1) = exp(sum(bi.*X(j,:),2))-prod((exp(bi)-1)./bi);
    
57.     end
    
58.     Y = Yn; %model outputs
    
59.     for j = 1:Nvar
    
60.        Xs(:,j) = sqrt(2)*erfinv(2.*X(:,j)-1); %model inputs
    
61.     end
    
62. elseif testfunc == 3 %Ishigami test case (non-linear/non-monotonic)
    
63.     Nvar = 3;
    
64.     lb1 = repmat([-pi -pi -pi],[Nevalpts 1]);
    
65.     ub1 = repmat([pi pi pi],[Nevalpts 1]);
    
66.     %Parameters
    
67.     a = 7; b = 0.1;
    
68.     f = @(x) (sin(x(:,1))+ a.*(sin(x(:,2)).^2)+(b.*(x(:,3).^4).*sin(x(:,1))));
    
69.     p = sobolset(Nvar,'skip',Delay);
    
70.     X = p(1:Nevalpts,:);
    
71.     Xss = lb1+((ub1-lb1).*X);
    
72.     %Evaluate model
    
73.     [Yn] = f(Xss);
    
74.     Y = Yn; %model output
    
75.     %Std normal space for Xs
    
76.     for j = 1:Nvar
    
77.         Xs(:,j) = sqrt(2)*erfinv(2.*X(:,j)-1); %model inputs
    
78.     end
    
79. end
    
80. 
    
81. %%Estimate alpha values for the Rolling Pin method
    
82. %using the maximum likelihood (MLE) optimization approach
    
83. %Here, an external call to LSHADEEpSinNLS algorithm is performed
    
84. %Note, on 2 cores this will likely take a bit of time, especially if
    
85. %Nevalpts is high
    
86. ry = Y;
    
87. rx = Xs;
    
88. Ngens = 20000; %This is the number of generations used for MLE optimization
    
89. Npop = 50; %This is the number of population members used for MLE optimization
    
90. nTol = 1E-6; %This is the tolerance used to moniter convergence
    
91. paramL = zeros(1,Nvar+1); paramU = ones(1,Nvar+1); %boundaries of the alpha parameter values, Rolling Pin method
    
92. paramU(1,1) = 0;
    
93. x2{1,1} = ry; x2{1,2} = rx; x2{1,3} = model;
    
94. [solset] = LSHADEpSinNLS(@MLEObjFuncFinalAll, Nvar+1, paramL, paramU, Npop, nTol, Ngens,x2); %Run optimization call
    
95. [m,p] = find(nonzeros(solset.Xbest(:,2)));
    
96. alpha = solset.Xbest(m(end),:); %Extract best alpha values for Rolling Pin method
    
97. 
    
98. %%Semi-Parametric method begins here
    
99. %Estimate Yi, transformed values, Equation 4
    
100. for j = 1:Nvar
     
101.     Yi(:,j) = (1-alpha(1,j+1)).*(ry(:,1))+alpha(1,j+1).*(rx(:,j));
     
102. end
     
103. 
     
104. %Find the PDFs and CDFs of the transformed distributions, Yi and Xs (Eqs.
     
105. %5 and 6)
     
106. parfor j = 1:Nvar
     
107.    tic
     
108.    pd1{j,1} = fitdist(Yi(:,j),'Kernel');
     
109.    FYr(:,j) = cdf(pd1{j,1},Yi(:,j));
     
110.    FX(:,j) = normcdf(rx(:,j),0,1);
     
111.    toc
     
112. end
     
113. 
     
114. %Fit Elliptical or Archimidean parametric copula model here using MATLAB's built in
     
115. %functions
     
116. for j = 1:Nvar
     
117.    u{j,1} = [FYr(:,j) FX(:,j)];
     
118.    %Construct Parametric Copula from CDF matrix, t-distribution is a
     
119.    %special case
     
120.    if model == 5
     
121.        [r{j,1},nus{j,1}] = copulafit(CopulaModel,u{j,1});
     
122.    else
     
123.        [r{j,1}] = copulafit(CopulaModel,u{j,1});
     
124.    end
     
125. end
     
126. 
     
127. %Sample from constructed parametric copula model, using MATLAB's built in
     
128. %functions
     
129. parfor j = 1:Nvar
     
130.     tic
     
131.     if model == 5
     
132.         cs = copularnd(CopulaModel,r{j,1},nus{j,1},Nevalpts);
     
133.     else
     
134.         cs = copularnd(CopulaModel,r{j,1},Nevalpts);
     
135.     end
     
136.     FYrn(:,j) = cs(:,1);
     
137.     FXn(:,j) = cs(:,2);
     
138.     Yin(:,j) = icdf(pd1{j,1},FYrn(:,j));
     
139.     Xn(:,j) = norminv(FXn(:,j),0,1);
     
140.     toc
     
141. end
     
142. 
     
143. %Find the marginal PDFs of the sampled dists, Yin and Xn (Eqs. 5 and 6)
     
144. parfor j = 1:Nvar
     
145.    tic
     
146.    FYpdf1(:,j) = pdf(pd1{j,1},Yin(:,j));
     
147.    Fpdf(:,j) = normpdf(Xn(:,j),0,1);
     
148.    toc
     
149. end
     
150. 
     
151. %Transform the sample back to original space, inverse of Equation 4
     
152. parfor j = 1:Nvar
     
153.    Ys(:,j) = (Yin(:,j)-alpha(1,j+1).*Xn(:,j))./(1-alpha(1,j+1));
     
154.    pd3{j,1} = fitdist(Ys(:,j),'Kernel');
     
155.    FyPDF(:,j) = pdf(pd3{j,1},Ys(:,j));
     
156. end
     
157. 
     
158. %Evaluate the parametric copula density
     
159. for j = 1:Nvar
     
160.     %Evaluate joint pdf directly at pts
     
161.     a = FYpdf1(:,j).*(1-alpha(1,j+1));
     
162.     b = normpdf(Xn(:,j),0,1).*(1-alpha(1,1));
     
163.     if model == 5
     
164.         fXa1(:,j) = copulapdf(CopulaModel,[FYrn(:,j),FXn(:,j)],r{j,1},nus{j,1});
     
165.     else
     
166.         fXa1(:,j) = copulapdf(CopulaModel,[FYrn(:,j),FXn(:,j)],r{j,1});
     
167.     end
     
168.     fX(:,j) = fXa1(:,j).*a.*b; %Equation 7
     
169. end
     
170. 
     
171. %Estimate the MI sensitivity indices, Equation 8
     
172. for j = 1:Nvar
     
173.     deltast(:,j) = abs(((FyPDF(:,j).*Fpdf(:,j))./(fX(:,j)))-1);
     
174.     delta(1,j) = 0.5*mean(deltast(:,j));
     
175. end
     
176. 
     
177. %Save results here
     
178. filename = strcat(['SemiParamCopula_',CopulaModel,'_Test',num2str(testfunc),'_N',num2str(Nevalpts),'_Rep',num2str(Rep),'.mat']);
     
179. save(filename,'delta','deltast','solset','Xs','Y','ry','rx');

复制代码

Semi-Parametric-Copula-Approach-master.zip (3.39 MB)

2019-9-20 18:14:21 上传

MixedVineToolbox-master.zip (41.03 KB)

2019-9-20 18:14:20 上传

扫码加我 拉你入群

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

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

分享0 收藏12 回帖

谢谢分享！

非常感謝您的分享！

joycelee727 发表于 2020-7-30 19:17
非常感謝您的分享！

能有用就好，这是源代码的。 要是你们能还创新就好了，修改作者的代码。哈哈

最近急用，非常感谢您的分享

已经修改售价为 0

您好，能否加上代码解析，看不懂啊

xkt510587 发表于 2021-6-12 19:47
您好，能否加上代码解析，看不懂啊

不太明白你说的， 需要我加什么代码解析？

tulipsliu 发表于 2019-9-20 18:15
已经下载论坛很多书籍，大概6000论坛币已经剩下只有300左右吧，此贴不收金币。我太累，已经看不了那么多书。 ...

感谢分享！！！

leeyaya 发表于 2021-6-18 16:33
感谢分享！！！

客气了。