matlab数据的基本统计分析

第四讲 数据的基本统计分析

数据的基本统计分析

1．数据的描述性统计分析

通常在得到数据并对数据进行除错的预处理后，需要对数据进行描述性的统计分析。比如：对数据中变量的最小值、最大值、中位数、平均值、标准差、偏度、峰度以及正态性检验等进行分析。对于这些经常性遇到的重复过程，我们可以自己编写函数，将函数保存在MATLAB自动搜索文件夹下，然后就可以直接调用自己定义的函数了。对于上述描述性统计分析，我们可以在MATLAB命令窗口中输入：edit discription，然后在弹出的窗口中选择yes，就创建了一个文件名为discription的M文件。然后在弹出的空白文件中编写以下M函数： function D=discription(x)

%descriptive statistic analysis

%input:

%x is a matrix, and each colummn stands for a variable

%output:

%D:structure variable,denotes Minimium,Maximium,Mean,Median,

%Standard_deviation,Skewness,Kurtosis,and normal distribution test,respectively.

%notes:when the number of oberservations of the colummn variables less than 30,

%Lilliefors test is used for normal distribution test,and output D.LSTA denotes

%test statistic and D.LCV denote critical value under 5% significant level;

%otherwise, Jarque-Bera test is used, and output D.JBSTA denotes test statistic

%and D.JBCV denote critical value under 5% significant level.If test statistic is

%less than critical value,the null hypothesis (normal distribution) can not

%be rejected under 5% significant level.

D.Minimium=min(x);

D.Maximium=max(x);

D.Mean=mean(x);

D.Median=median(x);

D.Standard_deviation=std(x);

D.Skewness=skewness(x);

D.Kurtosis=kurtosis(x);

if size(x,1)<30

disp('small observations,turn to Lilliefors test for normal distribution')

for i=1:size(x,2)

[h(i),p(i),Lilliefors(i),LCV(i)]=lillietest(x(:,i),0.05);

end

D.LSTA=Lilliefors;

D.LCV=LCV;

else

for i=1:size(x,2)

[h(i),p(i),Jarque_Bera(i),JBCV(i)]=jbtest(x(:,i),0.05);

end

D.JBSTA=Jarque_Bera;

D.JBCV=JBCV;

end

注意在上面给出的函数例子中，我们使用了discription作为文件名，这与函数文件中第一行中的discription保持了一致。这样就可以以D=discription(x)形式调用该函数。如果使用不同于discription的文件名保存，比如：statistic，则调用该函数时，必须以D=statistic(x)形式调用。为避免调用时的麻烦，尽量使用相同的名称保存函数。在上面的函数discription中给出了正态分布检验的统计量与5%显著水平下的临界值。当样本容量低于30时，使用Lilliefors 检验；当样本容量超过30时使用Jarque-Bera检验。

下面我们以上证综合指数为例来调用刚刚自定义的函数discription。假定我们只关心以开盘价、最高价、最低价、收盘价表示的日收益率。在读入数据并对数据进行除错的预处理后（将数据按照日期升序进行重新排列），我们得到变量b、c、d、e分别表示1990年12月19日到2006年9月27日之间的开盘价、最高价、最低价、收盘价数据。然后在MATLAB命令窗口中输入：

x=price2ret([b,c,d,e]);%将价格转换为对数收益率

D=discription(x)%调用自定义函数discription

得到以下结果：

D =

Minimium: [-0.3170 -0.1565 -0.4498 -0.1791]

Maximium: [0.7138 0.7607 0.7372 0.7192]

Mean: [7.4406e-004 7.3581e-004 7.4450e-004 7.3574e-004]

Median: [7.0916e-004 8.0367e-004 3.6515e-004 4.3624e-004]

Standard_deviation: [0.0291 0.0253 0.0278 0.0265]

Skewness: [4.5113 8.2876 4.2696 6.1913]

Kurtosis: [111.7483 229.2601 162.1498 156.0935]

JBSTA: [1.9186e+006 8.2927e+006 4.0928e+006 3.8010e+006]

JBCV: [5.9915 5.9915 5.9915 5.9915]

2．样本分布函数与概率密度函数

在对数据进行基本的描述性统计分析后，有时我们还需要对变量的样本分布函数与样本概率密度函数进行分析。甚至有时候，基于研究的需要，我们还要根据样本的历史数据，来产生随机样本进行某些研究。下面以1990年12月19日到2006年9月27日之间的上证综合指数收盘价为例，给出如何利用MATLAB得到上证综合指数日对数收益率的经验分布函数以及样本的概率密度函数，还有如何根据历史收益率的经验分布来生成随机数。

（1）样本分布函数

假定我们在MATLAB中已经读入了2000年1月1日到2006年6月1日之间的上证综合指数的日期和收盘价数据，在经过数据的预处理后，得到列向量a和

e，分别表示时期和收盘价。在MATLAB命令窗口下输入：

log_ret=price2ret(e);

h=figure;

set(h,'color','w')

plot(a(2:end),log_ret)

datetick('x',23)

xlabel('date')

ylabel('return')

title('daily return of Shanghai Composite')

图形输出结果如图所示。

上证综合指数日对数收益率

为了得到样本的分布函数，我们可以编写以下M函数，并以empirical_dist 的文件名保存在MATLAB自动搜索的文件夹下。

function [x,cumpr]=empirical_dist(data)

% generate empirical distribution function

% input:

% data is a vector

% output:

% x is sample observation vector

% cumpr is cumulative probability vector

if min(size(data))~=1

error('data must be a vector')

end

n=length(data);

data=reshape(data,n,1);

data=sort(data);

[x,a,b]=unique(data);

frequency=[a(1);diff(a)];

cumpr=cumsum(frequency)/n;

然后在MATLAB命令窗口下输入：

[x,cumpr]=empirical_dist(log_ret);

h=figure;

set(h,'color','w')

plot(x,cumpr)

ylabel('cumulative probability')

title('empirical distribution of daily returns on Shanghai Composite') 图形输出结果如图所示。

上证综合指数日对数收益率的经验分布

（2）样本概率密度函数

为了得到样本的概率密度函数，我们可以编写以下M函数，并以empirical_density的文件名保存在MATLAB自动搜索的文件夹下。

function [x,density]=empirical_density(data,m)

%generate relative frequency and probability density

%input:

%data is a vector

%m is number of intervals

% output:

% x is a vector points of intervals

% density is probability density

if min(size(data))~=1

error('data must be a vector')

end

n=length(data);

data=reshape(data,n,1);

zeta=min(abs(data))/10;

min1=min(data)-zeta;%locate low ending point

max1=max(data)+zeta;%locate high ending point

x=linspace(min1,max1,m+1);%generate intervals

density=hist(data,x)./(n*(x(2)-x(1)));

在上面的程序中，区间数目的由m确定。

利用前面得到的上证综合指数的日对数收益率log_ret，在MATLAB命令窗口下输入：

[x,density]=empirical_density(log_ret,200);

h1=figure(1);

set(h1,'color','w')

bar(x,hist(log_ret,x)/length(log_ret));

title('relative frequency');

h2=figure(2);

set(h2,'color','w')

plot(x,density);

title('probability density');

图形输出结果分别如图所示。

上证综合指数日对数收益率的相对频率

上证综合指数日对数收益率的概率密度函数

利用函数empirical_density有时候得到的概率密度函数显得不太光滑，这时，我们可以编写以下函数empirical_density1得到较为光滑概率密度函数。 function [x,density]=empirical_density1(data,m)

%generate smoothed probability density function

%input:

%data is a vector

%m is number of intervals

% output:

% x is a vector points of intervals

% density is smoothed probability density

if min(size(data))~=1

error('data must be a vector')

end

n=length(data);

data=reshape(data,n,1);

zeta=min(abs(data))/10;

min1=min(data)-zeta;%locate low ending point

max1=max(data)+zeta;%locate high ending point

x=linspace(min1,max1,m+1);%generate intervals

h=2*iqr(data)*length(data)^(-1/3);

density=zeros(1,length(x));

for j=1:n

density=density+normpdf(x,data(j),h);

end

density=density/n;

（3）产生服从经验分布的随机数

根据历史数据的经验分布产生随机数的基本原理如下：记()

f x为随机变量

F x服从[0，1]区间的均匀分布。X的概率密度函数，()

F x为其分布函数。则()

因此，先利用计算机产生服从[0，1]区间均匀分布的随机数y ，然后利用1()z F y ?=就得到了服从概率密度函数为()f x 的随机数z 。因此，如果知道()f x 或()F x ，应用上述原理就可以产生相应的随机数。

基于以上原理，我们可以编写以下函数名为generate_rv 的函数：

function random_number=generate_rv(data,m,n)

% generate random numbers based on empirical distribution

% input:

% data is a vector of sample points

% m is the number of rows of the matrix with generated random numbers (i.e.random_number)

% n is the number of columns of the matrix with generated random numbers % output:

% random_number is a m-by-n matrix with generated random numbers

[x,cumpr]=empirical_dist(data);

union_number=rand(m,n);

random_number=interp1(cumpr,x,union_number);

函数generate_rv 中用到了我们前面编写的函数empirical_dist。利用前面得到的上证综合指数的日对数收益率log_ret，在MATLAB 命令窗口下输入： h=figure;

set(h,'color','w')

random_number=generate_rv(log_ret,100,100);

plot(random_number,'.')

title('generated random numbers based on empirical distribution')

上述结果将根据上证综合指数日对数收益率的经验分布产生100行100列的相互独立的随机数。图形输出结果如图所示。

基于上证综合指数日对数收益率经验分布的随机数