In matlab or R, I want to solve the weight coefficient of multidimensional (three-dimensional or more) mixed copula, (not vine copula model).The mixture copula are constructed by clayton copula,gumbel copula and frank copula(or other copula), the following matlab EM algorithm code is only able to solve the weight coefficient of two-dimensional copula. Can someone teach me? Thank you very much!
clc;
data = xlsread('data_3-dimension-100set.xls');% The data can be randomly generated in
% the interval [0.1].
U=data(:,1);
V=data(:,2);
%The EM algorithm is used for estimation
S=3; % S is the number of Copula
N=length(U); % N is the sample size
th=0.000001; % Convergence condition
a=2.21;
gama=0.65; % Parameter initialization
rana=[1/3,1/3,1/3]'; % The weight parameter is initialized
a1=copulafit('Clayton',[U V]);
a2=copulafit('Gumbel',[U,V]);
a3=copulafit('Frank',[U,V]);
theta=[a1,a2,a3]'; % Initialization of dependent structure parameters
t=inf;
count=0; % Number of iterations
COPULA=zeros(N,S); %COPULA function matrix
ranack=zeros(S,1); % λ(k)*copula(n,k)
RZNK=zeros(N,S); %λ(k)*copula(n,k)/sum（λ(k)*copula(n,k)）
SCADD=zeros(S,1); %λ(k)*SCAD(λ(k))
turntheta_old=zeros(S,1);
linss1=zeros(3,51); %Used to store changes in the individual Copula weight parameters
linss2=zeros(3,51); %Used to store the variation of each Copula dependent parameter
while t>=th&&count<=50
rana_old=rana;
theta_old=theta;
turntheta_old(1)=log(theta_old(1));
turntheta_old(2)=log(theta_old(2)-1);
turntheta_old(3)=tan(theta_old(3));
%Since the dependent structure parameters are constrained and limited, the parameters are transformed so that unconstrained optimization methods are used to solve the
%Calculate the COPULA matrix
for n=1:N
COPULA(n,1)=copulapdf('Clayton',U(n,:),V(n,:),theta_old(1));
%Clayton Copula Density function value
COPULA(n,2)=copulapdf('Gumbel',U(n,:), V(n,:),theta_old(2));
%Gumble Copula Density function value
COPULA(n,3)=copulapdf('Frank',U(n,:), V(n,:),theta_old(3));
%Frank Copula Density function value
end
%computing SCADD
for k=1:S
SCADD(k)=rana_old(k)*SCAD(rana_old(k),a,gama);
end
delta=N*(1-sum(SCADD));
%E step
%rznk
for cn=1:N
for cm=1:S
ranack(cm)=rana_old(cm)*COPULA(cn,cm);
end
for cm=1:S
RZNK(cn,cm)=ranack(cm)/sum(ranack);
end
end
%M step
% Find the update equation for the weight parameters
for cm=1:S
rana(cm)=(sum(RZNK(:,cm))-N*rana_old(cm)*SCAD(rana_old(cm),a,gama))/delta;
end
if sum(rana<0)>=1
rana=rana_old;
break;%Determine if there is a situation where rana is less than 0, then terminate the jump out loop
end
rana;
linss1(:,count+1)=rana;
%Finding update equations for dependent structural parameters
%Updating equations for dependent structural parameters using the BFGS approach
[turntheta,fval,exitflag,output,grad,hessian]=runtheta(N,U,V,RZNK,turntheta_old);
turntheta;
%Convert the dependent structural parameters back
theta(1)=exp(turntheta(1));
theta(2)=exp(turntheta(2))+1;
%Preventing the emergence of Frank Copula's dependence structure parameterized by 0
if turntheta(3)==0
theta(3)=turntheta(3)+0.0001;
else
theta(3)=atan(turntheta(3));
end
theta;
linss2(:,count+1)=theta;
%Judgment termination conditions
t=max([norm(rana_old(:)-rana(:));norm(theta_old(:)-theta(:))]);
count=count+1;
end
disp(rana)
disp(theta)
%%% runtheta.m
function [turntheta,fval,exitflag,output,grad,hessian]=runtheta(N,U,V,RZNK,theta)
options=optimset('TolFun',1e-
10,'LargeScale','off','FinDiffType','central','HessUpdate','bfgs');
[turntheta,fval,exitflag,output,grad,hessian]=fminunc(@mtheta,theta,options);
function y=mtheta(theta)
%theta :dependent structural parameters
%U,V: data
%rana_old: Structural parameters of theta
%rznk为入(k)*copula(n，k)/sum(入(k)*copula(n，k))
%N : sample
if theta(3)==0
theta(3)=theta(3)+0.0001;
end
y=0;
for t=1:N
y=y+(log((1+(exp(theta(1))))*(U(t,1)*V(t,1))^(-(exp(theta(1)))-1)*(U(t,1)^(-(exp(theta(1))))+V(t,1)^(-(exp(theta(1))))-1)^(-2-1/(exp(theta(1)))))*RZNK(t,1)+...
log(exp(-((-log(U(t,1))).^(exp(theta(2))+1)+(-log(V(t,1))).^(exp(theta(2))+1)).^(1/(exp(theta(2))+1)))*(((-log(U(t,1))).^(exp(theta(2))+1)+(-log(V(t,1))).^(exp(theta(2))+1)).^(1/(exp(theta(2))+1))+(exp(theta(2))+1)-1)*(log(U(t,1))*log(V(t,1)))^((exp(theta(2))+1)-1)/(U(t,1)*V(t,1)*((-log(U(t,1))).^(exp(theta(2))+1)+(-log(V(t,1))).^(exp(theta(2))+1))^(2-1/(exp(theta(2))+1))))*RZNK(t,2)+...
log(theta(3)*(1-exp(-theta(3)))*exp(-theta(3)*(U(t,1)+V(t,1)))/((1-exp(-theta(3)))-(1-exp(-theta(3)*U(t,1)))*(1-exp(-theta(3)*V(t,1))))^2)*RZNK(t,3));
end
y=-y;
end
end
%%%SCAD.m
function SCAD=SCAD(rana,a,gama) %derivative of the SCAD function
if rana<0||a<=2
display('Parameter error!')
return;
end
if rana==0
SCAD=0;
elseif gama-rana>=0
SCAD=gama;
else
if a*gama-rana>0
SCAD=(a*gama-rana)/(a-1);
else
SCAD=0;
end
end
end