Copula理论及MATLAB应用实例.docx
《Copula理论及MATLAB应用实例.docx》由会员分享，可在线阅读，更多相关《Copula理论及MATLAB应用实例.docx（35页珍藏版）》请在冰豆网上搜索。
Copula理论及MATLAB应用实例
%--------------------------------------------------------------------------
%Copula理论及其在matlab中的实现程序应用实例
%--------------------------------------------------------------------------
%******************************读取数据*************************************
%从文件hushi.xls中读取数据
hushi=xlsread（'hushi.xls'）;
%提取矩阵hushi的第5列数据，即沪市的日收益率数据
X=hushi（:
5）;
%从文件shenshi.xls中读取数据
shenshi=xlsread（'shenshi.xls'）;
%提取矩阵shenshi的第5列数据，即深市的日收益率数据
Y=shenshi（:
5）;
%****************************绘制频率直方图*********************************
%调用ecdf函数和ecdfhist函数绘制沪、深两市日收益率的频率直方图
[fx,xc]=ecdf（X）;
figure;
ecdfhist（fx,xc,30）;
xlabel（'沪市日收益率'）;%为X轴加标签
ylabel（'f（x）'）;%为Y轴加标签
[fy,yc]=ecdf（Y）;
figure;
ecdfhist（fy,yc,30）;
xlabel（'深市日收益率'）;%为X轴加标签
ylabel（'f（y）'）;%为Y轴加标签
%****************************计算偏度和峰度*********************************
%计算X和Y的偏度
xs=skewness（X）
ys=skewness（Y）
%计算X和Y的峰度
kx=kurtosis（X）
ky=kurtosis（Y）
%******************************正态性检验***********************************
%分别调用jbtest、kstest和lillietest函数对X进行正态性检验
[h,p]=jbtest（X）%Jarque-Bera检验
[h,p]=kstest（X,[X,normcdf（X,mean（X）,std（X））]）%Kolmogorov-Smirnov检验
[h,p]=lillietest（X）%Lilliefors检验
%分别调用jbtest、kstest和lillietest函数对Y进行正态性检验
[h,p]=jbtest（Y）%Jarque-Bera检验
[h,p]=kstest（Y,[Y,normcdf（Y,mean（Y）,std（Y））]）%Kolmogorov-Smirnov检验
[h,p]=lillietest（Y）%Lilliefors检验
%****************************求经验分布函数值*******************************
%调用ecdf函数求X和Y的经验分布函数
[fx,Xsort]=ecdf（X）;
[fy,Ysort]=ecdf（Y）;
%调用spline函数，利用样条插值法求原始样本点处的经验分布函数值
U1=spline（Xsort（2:
end）,fx（2:
end）,X）;
V1=spline（Ysort（2:
end）,fy（2:
end）,Y）;
%调用ecdf函数求X和Y的经验分布函数
[fx,Xsort]=ecdf（X）;
[fy,Ysort]=ecdf（Y）;
%提取fx和fy的第2个至最后一个元素，即排序后样本点处的经验分布函数值
fx=fx（2:
end）;
fy=fy（2:
end）;
%通过排序和反排序恢复原始样本点处的经验分布函数值U1和V1
[Xsort,id]=sort（X）;
[idsort,id]=sort（id）;
U1=fx（id）;
[Ysort,id]=sort（Y）;
[idsort,id]=sort（id）;
V1=fy（id）;
%*******************************核分布估计**********************************
%调用ksdensity函数分别计算原始样本X和Y处的核分布估计值
U2=ksdensity（X,X,'function','cdf'）;
V2=ksdensity（Y,Y,'function','cdf'）;
%**********************绘制经验分布函数图和核分布估计图**********************
[Xsort,id]=sort（X）;%为了作图的需要，对X进行排序
figure;%新建一个图形窗口
plot（Xsort,U1（id）,'c','LineWidth',5）;%绘制沪市日收益率的经验分布函数图
holdon
plot（Xsort,U2（id）,'k-.','LineWidth',2）;%绘制沪市日收益率的核分布估计图
legend（'经验分布函数','核分布估计','Location','NorthWest'）;%加标注框
xlabel（'沪市日收益率'）;%为X轴加标签
ylabel（'F（x）'）;%为Y轴加标签
[Ysort,id]=sort（Y）;%为了作图的需要，对Y进行排序
figure;%新建一个图形窗口
plot（Ysort,V1（id）,'c','LineWidth',5）;%绘制深市日收益率的经验分布函数图
holdon
plot（Ysort,V2（id）,'k-.','LineWidth',2）;%绘制深市日收益率的核分布估计图
legend（'经验分布函数','核分布估计','Location','NorthWest'）;%加标注框
xlabel（'深市日收益率'）;%为X轴加标签
ylabel（'F（x）'）;%为Y轴加标签
%****************************绘制二元频数直方图*****************************
%调用ksdensity函数分别计算原始样本X和Y处的核分布估计值
U=ksdensity（X,X,'function','cdf'）;
V=ksdensity（Y,Y,'function','cdf'）;
figure;%新建一个图形窗口
%绘制边缘分布的二元频数直方图，
hist3（[U（:
）V（:
）],[30,30]）
xlabel（'U（沪市）'）;%为X轴加标签
ylabel（'V（深市）'）;%为Y轴加标签
zlabel（'频数'）;%为z轴加标签
%****************************绘制二元频率直方图*****************************
figure;%新建一个图形窗口
%绘制边缘分布的二元频数直方图，
hist3（[U（:
）V（:
）],[30,30]）
h=get（gca,'Children'）;%获取频数直方图的句柄值
cuv=get（h,'ZData'）;%获取频数直方图的Z轴坐标
set（h,'ZData',cuv*30*30/length（X））;%对频数直方图的Z轴坐标作变换
xlabel（'U（沪市）'）;%为X轴加标签
ylabel（'V（深市）'）;%为Y轴加标签
zlabel（'c（u,v）'）;%为z轴加标签
%***********************求Copula中参数的估计值******************************
%调用copulafit函数估计二元正态Copula中的线性相关参数
rho_norm=copulafit（'Gaussian',[U（:
）,V（:
）]）
%调用copulafit函数估计二元t-Copula中的线性相关参数和自由度
[rho_t,nuhat,nuci]=copulafit（'t',[U（:
）,V（:
）]）
%********************绘制Copula的密度函数和分布函数图************************
[Udata,Vdata]=meshgrid（linspace（0,1,31））;%为绘图需要，产生新的网格数据
%调用copulapdf函数计算网格点上的二元正态Copula密度函数值
Cpdf_norm=copulapdf（'Gaussian',[Udata（:
）,Vdata（:
）],rho_norm）;
%调用copulacdf函数计算网格点上的二元正态Copula分布函数值
Ccdf_norm=copulacdf（'Gaussian',[Udata（:
）,Vdata（:
）],rho_norm）;
%调用copulapdf函数计算网格点上的二元t-Copula密度函数值
Cpdf_t=copulapdf（'t',[Udata（:
）,Vdata（:
）],rho_t,nuhat）;
%调用copulacdf函数计算网格点上的二元t-Copula分布函数值
Ccdf_t=copulacdf（'t',[Udata（:
）,Vdata（:
）],rho_t,nuhat）;
%绘制二元正态Copula的密度函数和分布函数图
figure;%新建图形窗口
surf（Udata,Vdata,reshape（Cpdf_norm,size（Udata）））;%绘制二元正态Copula密度函数图
xlabel（'U'）;%为X轴加标签
ylabel（'V'）;%为Y轴加标签
zlabel（'c（u,v）'）;%为z轴加标签
figure;%新建图形窗口
surf（Udata,Vdata,reshape（Ccdf_norm,size（Udata）））;%绘制二元正态Copula分布函数图
xlabel（'U'）;%为X轴加标签
ylabel（'V'）;%为Y轴加标签
zlabel（'C（u,v）'）;%为z轴加标签
%绘制二元t-Copula的密度函数和分布函数图
figure;%新建图形窗口
surf（Udata,Vdata,reshape（Cpdf_t,size（Udata）））;%绘制二元t-Copula密度函数图
xlabel（'U'）;%为X轴加标签
ylabel（'V'）;%为Y轴加标签
zlabel（'c（u,v）'）;%为z轴加标签
figure;%新建图形窗口
surf（Udata,Vdata,reshape（Ccdf_t,size（Udata）））;%绘制二元t-Copula分布函数图
xlabel（'U'）;%为X轴加标签
ylabel（'V'）;%为Y轴加标签
zlabel（'C（u,v）'）;%为z轴加标签
%**************求Kendall秩相关系数和Spearman秩相关系数***********************
%调用copulastat函数求二元正态Copula对应的Kendall秩相关系数
Kendall_norm=copulastat（'Gaussian',rho_norm）
%调用copulastat函数求二元正态Copula对应的Spearman秩相关系数
Spearman_norm=copulastat（'Gaussian',rho_norm,'type','Spearman'）
%调用copulastat函数求二元t-Copula对应的Kendall秩相关系数
Kendall_t=copulastat（'t',rho_t）
%调用copulastat函数求二元t-Copula对应的Spearman秩相关系数
Spearman_t=copulastat（'t',rho_t,'type','Spearman'）
%直接根据沪、深两市日收益率的原始观测数据，调用corr函数求Kendall秩相关系数
Kendall=corr（[X,Y],'type','Kendall'）
%直接根据沪、深两市日收益率的原始观测数据，调用corr函数求Spearman秩相关系数
Spearman=corr（[X,Y],'type','Spearman'）
%******************************模型评价*************************************
%调用ecdf函数求X和Y的经验分布函数
[fx,Xsort]=ecdf（X）;
[fy,Ysort]=ecdf（Y）;
%调用spline函数，利用样条插值法求原始样本点处的经验分布函数值
U=spline（Xsort（2:
end）,fx（2:
end）,X）;
V=spline（Ysort（2:
end）,fy（2:
end）,Y）;
%定义经验Copula函数C（u,v）
C=@（u,v）mean（（U<=u）.*（V<=v））;
%为作图的需要，产生新的网格数据
[Udata,Vdata]=meshgrid（linspace（0,1,31））;
%通过循环计算经验Copula函数在新产生的网格点处的函数值
fori=1:
numel（Udata）
CopulaEmpirical（i）=C（Udata（i）,Vdata（i））;
end
figure;%新建图形窗口
%绘制经验Copula分布函数图像
surf（Udata,Vdata,reshape（CopulaEmpirical,size（Udata）））
xlabel（'U'）;%为X轴加标签
ylabel（'V'）;%为Y轴加标签
zlabel（'EmpiricalCopulaC（u,v）'）;%为z轴加标签
%通过循环计算经验Copula函数在原始样本点处的函数值
CUV=zeros（size（U（:
）））;
fori=1:
numel（U）
CUV（i）=C（U（i）,V（i））;
end
%计算线性相关参数为0.9264的二元正态Copula函数在原始样本点处的函数值
rho_norm=0.9264;
Cgau=copulacdf（'Gaussian',[U（:
）,V（:
）],rho_norm）;
%计算线性相关参数为0.9325，自由度为4的二元t-Copula函数在原始样本点处的函数值
rho_t=0.9325;
k=4.0089;
Ct=copulacdf（'t',[U（:
）,V（:
）],rho_t,k）;
%计算平方欧氏距离
dgau2=（CUV-Cgau）'*（CUV-Cgau）
dt2=（CUV-Ct）'*（CUV-Ct）
灰色预测[GM（1,1）]MATLAB程序
%本程序主要用来计算根据灰色理论建立的模型的预测值。
%应用的数学模型是GM（1,1）。
%原始数据的处理方法是一次累加法。
y=input（'请输入数据'）;%输入数据请用如例所示形式：
[48.757.1768.7692.15]
n=length（y）;
yy=ones（n,1）;
yy
（1）=y
（1）;
fori=2:
n
yy（i）=yy（i-1）+y（i）;
end
B=ones（n-1,2）;
fori=1:
（n-1）
B（i,1）=-（yy（i）+yy（i+1））/2;
B（i,2）=1;
end
BT=B';
forj=1:
n-1
YN（j）=y（j+1）;
end
YN=YN';
A=inv（BT*B）*BT*YN;
a=A
（1）;
u=A
（2）;
t=u/a;
t_test=input（'请输入需要预测个数：
'）;
i=1:
t_test+n;
yys（i+1）=（y
（1）-t）.*exp（-a.*i）+t;
yys
（1）=y
（1）;
forj=n+t_test:
-1:
2
ys（j）=yys（j）-yys（j-1）;
end
x=1:
n;
xs=2:
n+t_test;
yn=ys（2:
n+t_test）;
plot（x,y,'^r',xs,yn,'*-b'）;
det=0;
fori=2:
n
det=det+abs（yn（i）-y（i））;
end
det=det/（n-1）;
disp（['百分绝对误差为：
',num2str（det）,'%']）;
disp（['预测值为：
',num2str（ys（n+1:
n+t_test））]）;
Matlab二维绘图
2010-06-0820:
41
本节介绍MATLAB的两种基本绘图功能：
二维平面图形和三维立体图形。
5.1二维平面图形
5.1.1基本图形函数
plot是绘制二维图形的最基本函数，它是针对向量或矩阵的列来绘制曲线的。
也就是
说，使用plot函数之前，必须首先定义好曲线上每一点的x及y坐标，常用格式为：
（1）plot（x）当x为一向量时，以x元素的值为纵坐标，x的序号为横坐标值绘制
曲线。
当x为一实矩阵时，则以其序号为横坐标，按列绘制每列元素值相对于其序号的曲
线，当x为m×n矩阵时，就由n条曲线。
（2）plot（x,y）以x元素为横坐标值，y元素为纵坐标值绘制曲线。
（3）plot（x,y1,x,y2,…）以公共的x元素为横坐标值，以y1,y2,…元素为纵坐标值绘制多条曲线。
例5.1.1画出一条正弦曲线和一条余弦曲线。
>>x=0:
pi/10:
2*pi;
>>y1=sin（x）;
>>y2=cos（x）;
>>plot（x,y1,x,y2）
图5.1.1函数plot绘制的正弦曲线
在绘制曲线图形时，常常采用多种颜色或线型来区分不同的数据组，MATLAB软件专
门提供了这方面的参数选项（见表5.1.1），我们只要在每个坐标后加上相关字符串，就可
实现它们的功能。
表5.1.1绘图参数表
色彩字符颜色线型字符线型格式标记符号数据点形式标记符号数据点形式
y黄-实线.点<小于号
m紫：
点线o圆s正方形
c青-.点划线x叉号d菱形
r红--虚线+加号h六角星
g绿*星号p五角星
b蓝v向下三角形
w白^向上三角形
k黑>大于号
例如，在上例中输入
>>plot（x,y1,'r+-',x,y2,'k*:
'）
则得图5.1.2
图5.1.2使用不同标记的plot函数绘制的正弦曲线
5.1.2图形修饰
MATLAB软件为用户提供了一些特殊的图形函数，用于修饰已经绘制好的图形。
表5.1.2图形修饰函数表
函数 含义
gridon（/off） 给当前图形标记添加（取消）网络
xlable（‘string’） 标记横坐标
ylabel（‘string’） 标记纵坐标
title（‘string’） 给图形添加标题
text（x,y,’string’） 在图形的任意位置增加说明性文本信息
gtext（‘string’） 利用鼠标添加说明性文本信息
axis（[xminxmaxyminymax]）设置坐标轴的最小最大值
例5.1.2给例5.1.1的图形中加入网络和标记。
（见图5.1.3和5.1.4）
>>x=0:
pi/10:
2*pi;
>>y1=sin（x）;
>>y2=cos（x）;
>>plot（x,y1,x,y2）
>>gridon
>>xlabel（'independentvariableX'）
>>ylabel（'DependentVariableY1&Y2'）
>>title（'SineandCosineCurve'）
>>text（1.5,0.3,'cos（x）'）
>>gtext（'sin（x）'）
>>axis（[02*pi-0.90.9]）
图5.1.3使用了图形修饰的plot函数绘制的正弦曲线
5.1.3图形的比较显示
在一般默认的情况下，MATLAB每次使用plot函数进行图形绘制，将重新产生一个图
形窗口。
但有时希望后续的图形能够和前面所绘制的图形进行比较。
一般来说有两种方法：
一是采用holdon（/off）命令，将新产生的图形曲线叠加到已有的图形上；
二是采用subplot（m,n,k）函数，将图形窗口分隔成nm×个子图，并选择第k个子图作为当前图形
，然后在同一个视图窗口中画出多个小图形。
例5.1.3在同一窗口中绘制线段。
>>x=0:
pi/10:
2*pi;
>>y1=sin（x）;
>>y2=cos（x）;
>>y3=x;
>>y4=log（x）;
>>plot（x,y1,x,y2）
>>holdon
>>plot（x,y3）
>>plot（x,y4）
>>holdoff
例5.1.4在多个窗口中绘制图形。
（见图5.1.6）
>>x=0:
pi/10:
2*pi;
>>y1=sin（x）;
>>y2=cos（x）;
>>y3=exp（x）;
>>y4=log（x）;
>>subplot（2,2,1）;
>>plot（x,y1）;
>>subplot（2,2,2）;
>>plot（x,y2）;
>>subplot（2,2,3）;
>>plot（x,y3）;
>>subplot（2,2,4）;
>>plot（x,y4）;
[说明]
（1）子窗口的序号按行由上往下，按列从左向右编号。
（2）如果不用指令clf清除，以后图形将被绘制在子图形窗口中。
图5.1.6图形的比较显示（图形窗口分割方法）
5.2三维立体图形
5.2.1三维曲线图
与二维图形相对应，MATLAB提供了plot3函数，可以在三维空间中绘制三维曲线，
它的格式类似于plot，不过多了z方向的数据。
plot3的调用格式为：
plot3（x1,y1,z1,x2,y2,z2,...）
其中x1,y1,z1,x2,y2,z2,…等分别为维数相同的向量，分别存储着曲线的三个坐标值，该
函数的
升级会员 