极大似然估计

金融计算与编程 上海财经大学金融学院 曹志广 1 极大似然估计极大似然估计极大似然估计极大似然估计 极大似然估计方法在金融领域中的应用十分广泛。该方法利用已知的概率密 度函数形式，构造对数似然函数，然后最大化该似然函数从而求得概率密度函数 中所含的参数估计量。比如：对 GARCH(1,1)模型中的参数估计中，如果均值方 程中的扰动项服从正态分布，则我们可以利用正态分布的概率密度函数对所含参 数进行估计。 1．极大似然估计基本原理 (1)参数估计 下面以线性回归中系数的极大似然估计为例来说明极大似然估计基本原理。 考虑线性回归：Y X β ε= + ， 2~ (0, )Y X Nε β σ= − 则对于 X 和Y 的每一对观测值 ( , )i iX Y ，这里， iX 为行向量，其概率密度函 数形式如下： 2 2 1 1( , ) exp( ( ) ) 22 i i i i Y Xf X Y β σpiσ − = − 给定 N 对相互独立的观测值 ( , )i iX Y ， 1,2,...,i N= ，样本中所有观测值的总体概 率密度函数 ( , )L β σ 为单个观测值概率密度函数的乘积，即： 2 2 1 1 1( , ) exp( ( ) ) 22 N i i i Y XL ββ σ σpiσ= − = −∏ （1） 极大似然估计要给出参数 ( , )β σ 的估计量使得（1）式最大。由于（1）式为乘积 的形式，直接对最大化（1）式求解最优解，比较麻烦。因此，采用似然函数的 对数形式： 2 22 1 1 1( , ) [ ( ) ( ) ] 22 N i i i LnL Ln Y Xβ σ β σpiσ= = − −∑ 然后求解以下最优化问 题 快递公司问题件快递公司问题件货款处理关于圆的周长面积重点题型关于解方程组的题及答案关于南海问题 ： 2 22( , ) 1 1 1 max ( , ) [ ( ) ( ) ] 22 N i i i LnL Ln Y X β σ β σ β σpiσ= = − −∑ （2） 最后得到的参数 ( , )β σ 的估计量与普通最小二乘法得到的结果一样。因此，当普 通最小二乘法回归方程中的残差服从正态分布时，普通最小二乘估计与极大似然 估计的结果是一样的。 更一般地，我们用θ 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 示需要估计的参数向量，相应地对数似然函数为： ( )LnL θ 。 (2)参数估计的 标准 excel标准偏差excel标准偏差函数exl标准差函数国标检验抽样标准表免费下载红头文件格式标准下载 误差 求解优化问题（2），虽然给出了参数θ 的估计量 ˆθ ，但并没有给出估计的标 金融计算与编程 上海财经大学金融学院 曹志广 2 准误差。如果对数似然函数 ( )LnL θ 在其估计量 ˆθ 处的二阶倒数的期望是已知的， 则极大似然估计量的渐进协方差矩阵 1[ ( )]I θ − 满足： 2 1 1 1( ) ( ) ( )[ ( )] { [ ]} { [( )( )]}LnL LnL LnLI E Eθ θ θθ θ θ θ θ − − − ∂ ∂ ∂ = − = ′ ′∂ ∂ ∂ ∂ （3） 通常情况下 ( )LnL θ 是一个非常复杂的非线性函数，我们很难得到（3）式中期望 值的解析解形式。因此，根据（3）式来求解极大似然估计量的渐进协方差矩阵 1[ ( )]I θ − 的解析解形式很困难。其中一个方法就是利用下面的式子（4）来估计 1[ ( )]I θ − ，其中 1ˆˆ[ ( )]I θ − 为 1[ ( )]I θ − 的估计。 1 1 1 ˆˆ ˆ ˆ[ ( )] [ ] N i i i I g gθ − − = ′= ∑ （4） 其中： ˆ[ ( , )] ˆ ˆ i i Ln f X g θ θ ∂ = ∂ [ ( , )]iLn f X θ 表示在观测值 ( , )i iX Y 处的对数似然函数。 由于（4）式中的 [ ( , )]iLn f X θ 在通常情况下仍然是一个复杂的非线性函数， 求解 ˆig 的解析形式也非常困难。因此，我们通常以数值解形式（5）来近似 ˆig 。 ˆ ˆ ˆ ˆ[ ( , )] [ ( , )] ˆ ˆ2 i i i i Ln f X Ln f X g h θ θ θ θ θ + ∆ − − ∆ ≈ = ∆ （5） 因此，极大似然估计量的渐进协方差矩阵估计量为 1 1 [ ] N i i i h h − = ′∑ 。 2．极大似然估计的 MATLAB 函数 基于以上极大似然估计的原理，我们可以编写以下 MATLAB 函数进行极大似 然估计。 function [para,standard_deviation,fv]=my_mle(fun,para0,varargin) %estimate parameters and standard errors when using maximium likelihood estimation(MLE) %input %fun: a function defined by users for calculating log probability density function (pdf) and negative sum of logarithm of pdf %para0 : given initial parameters % varargin: other needed parameters required by fun 金融计算与编程 上海财经大学金融学院 曹志广 3 %output %para: estimated parameters % standard_deviation: standard deviations of estimated parameters % fv: maximized likelihood function value %%%%%%%%%%% %example1: estimate mean and standard deviation by realizations of a %random variable which is normally distributed %function f=mynormpdfsum(x,num,y) %yy=1/sqrt(2*pi)/x(2)*exp(-(y-x(1)).^2/2/x(2)^2); %if num==1 %(note: it must be set to 1) %f=log(yy); %else f=-sum(log(yy));end %%%%%%%%% %y=2+3*randn(5000,1); %[para,standard_deviation]=my_mle('mynormpdfsum',[0;2],y) %%%%%%%%% % example2:estimate coefficients in a linear regression % clear % x=randn(500,1); % y=2+3*x+randn(500,1); % [para,standard_deviation,fv]=my_mle('mynormpdfsum001',[1;2;3],y,x) % %%%%%%%%%%% % function f=mynormpdfsum001(para,num,y,x) % yy=1/sqrt(2*pi)/para(3)*exp(-(y-para(1)-para(2)*x).^2/2/para(3)^2); % if num==1 %(note: it must be set to 1) % f=log(yy); % else f=-sum(log(yy));end %%%%%%%%%% %Zhiguang Cao, 2006,2,23 %caozhiguang@21cn.com para0=para0(:); [para,fv]=fminsearch(fun,para0,[],2,varargin{:}); fv=-fv; n=length(para); d=numericalfirstderivative(fun,para,1,varargin{:}); standard_deviation=sqrt(diag(pinv(d'*d))); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function f=numericalfirstderivative(fun,parameter,varargin) % input: % fun: the name of a function % parameter: given parameter with respect to which first-order derivative is calculated % varargin: other needed parameters required by fun 金融计算与编程 上海财经大学金融学院 曹志广 4 % output: % f: numerical first order derivative of fun at parameter n=length(parameter); for i=1:n a=zeros(n,1); a(i)=min(parameter(i)*1e-6,1e-5); y1(:,i)=feval(fun,parameter+a,varargin{:}); y2(:,i)=feval(fun,parameter-a,varargin{:}); f(:,i)=(y1(:,i)-y2(:,i))/2/a(i); end 下面我们给出一个例子。考虑一个随机变量 X ， X 以α 的概率取 1X ， 1X 服 从均值为 1µ ，方差为 21σ 的正态分布；以1 α− 的概率取 2X ， 2X 服从均值为 2µ ， 方差为 2 2σ 的正态分布。现假定观测到随机变量 X 的 2000 个实现值，试估计参数 1 2 1 2( , , , , )θ α µ µ σ σ ′= 。 我们得到单个观测值 ix 的概率密度函数为： 2 21 2 2 2 1 21 2 1 1 1 1( , ) exp( ( ) ) (1 ) exp( ( ) ) 2 22 2 i i i x xf x µ µθ α α σ σpiσ piσ − − = − + − − 我们首先编写以下形式的函数 mle_ex_001： function f=mle_ex_001(parameter,num,observations) alpha=parameter(1);mean1=parameter(2);mean2=parameter(3); std1=parameter(4);std2=parameter(5); y=alpha*1/sqrt(2*pi)/std1*exp(-(observations-mean1).^2/2/std1^2)+... (1-alpha)*1/sqrt(2*pi)/std2*exp(-(observations-mean2).^2/2/std2^2); if num==1 f=log(y); else f=-sum(log(y)); end 然后在 MATLAB主窗口下输入： clear n=2000;% set number of observations alpha=0.7;mean1=0;mean2=1;std1=1;std2=2;%set parameters for i=1:n rv=rand(1); %%%%%%%% %generate normally distributed random variable with mean 0 and standard %deviation 1 with probability 0.7. if rv<=alpha 金融计算与编程 上海财经大学金融学院 曹志广 5 observations(i)=randn(1); else %generate normally distributed random variable with mean 1 and standard %deviation 2 with probability 0.3. observations(i)=normrnd(1,2); end end observations=observations(:); [para,standard_deviation]=my_mle('mle_ex_001',[0.5;0;0;1;1],observations) 得到结果如下： para = 0.3038 1.2578 -0.0344 1.9639 0.9719 standard_deviation = 0.0446 0.1730 0.0421 0.0827 0.0420 即得到该随机变量以 30.38%的概率为服从均值为 1.2578，标准差为 1.9639的正 态分布；以 69.62%的概率为服从均值为-0.0344，标准差为 0.9719的正态分布。 这个结果与我们设定的以 30%的概率为服从均值为 1，标准差为 2 的正态分布； 以 70%的概率为服从均值为 0，标准差为 1 的正态分布比较接近。 3333．．．．上证综合指数收益率广义双曲线分布的极大似然估计上证综合指数收益率广义双曲线分布的极大似然估计上证综合指数收益率广义双曲线分布的极大似然估计上证综合指数收益率广义双曲线分布的极大似然估计 金融资产收益率的分布，对金融资产投资、风险管理等具有重要意义，吸引 了众多学者研究这个问题。现实金融数据的分布通常表现为厚尾性和不对称性， 因此用正态分布拟合实际金融数据的分布有很大的局限性，比如在 VAR 的计算 中，由于金融数据分布的厚尾性，在正态分布的假设条件下计算 VAR会带来较大 的误差。为此许多学者开始寻求更为合理的分布假设。Mandelbrot 提出了用稳 定分布代替金融数据正态分布的假设，但稳定分布的尾部通常比实际分布要更 厚。又有学者提出用截尾的稳定分布作为证券收益率的分布，截尾的稳定分布实 际上是中间部分仍用稳定分布, 两个尾部用比负幂律分布瘦的指数分布来代替 的一种混合分布，但运用截尾的稳定分布假设时，应当在何处截尾是一个问题。 因此，又有许多学者开始转向广义双曲线分布（Generalized Hyperbolic Distribution）。1977 年 Barndorff-Nielsen提出了广义双曲线分布，Eberlein 和 Keller(1995)则首先将其应用到了金融领域，由于广义双曲线分布的尾部要 比稳定分布的尾部要“薄”，因此广义双曲线分布在金融领域中得到了迅速发展。 广义双曲线分布的概率密度函数定义如下： )(22 5.0 25.02/22 ))(())()(,,,,(),,,,,( µβλλ µδαµδµδβαλµδβαλ −−− −+−+= xexKxgxGH 其中： (.)νK 为调整后的第二类 Bessel 函数。 金融计算与编程 上海财经大学金融学院 曹志广 6 2 2 / 2 ( 0.5) 2 2 ( )( , , , , ) 2 ( ) 0 0, ; 0 0, ; 0 0, g K λ λ λ λ α βλ α β δ µ piα δ δ α β λ δ β α λ δ β α λ δ β α − − = − > ≥ < = > < < > ≤当 时， 当 时， 当 时， 广义双曲线分布可以派生出以下两个分布：逆高斯分布（当 5.0−=λ ）和双 曲线分布（当 1=λ ）。 下面以 1990年 12月 20日—2005年 12月 31日的上证综指日对数收益率为 例进行了广义双曲线分布的参数估计。 假定我们在MATLAB主窗口下已经得到上证综指日对数收益率y，则在 MATLAB 主窗口下输入： [para,standard_deviation]=my_mle('gh_log_fun',[-1;25;1;0.0005;0.02],y) 得到结果如下： para = -0.8206 8.1582 0.7939 0.0001 0.0134 standard_deviation = 0.0674 2.0768 0.8736 0.0003 0.0008 函数 gh_log_fun 的具体形式如下： function f=gh_log_fun(x,num,y) lambda=x(1);alpha=x(2);beta=x(3);mu=x(4);delta=x(5); a=(alpha^2-beta^2)^(lambda/2); b=sqrt(2*pi)*alpha.^(lambda-0.5)*delta^(lambda); c=besselk(lambda,delta*sqrt(alpha.^2-beta^2)); g=a/(b*c); a1=(delta^2+(y-mu).^2).^(lambda/2-0.25); b1=besselk(lambda-0.5,alpha*sqrt(delta^2+(y-mu).^2)); c1=exp(beta*(y-mu)); y=g*a1.*b1.*c1; if num==1 f=log(y); else f=-sum(log(y)); end 金融计算与编程 上海财经大学金融学院 曹志广 7 然后在 MATLAB主窗口下输入： [x,density]=empirical_density(y,200); [xx,pdf]=gh_pdf(para,y); h=figure; set(h,'color','w') plot(xx,pdf,'r-') hold on plot(x,density,'b:') legend('generalized hyperbolic pdf','sample pdf',1) 结果得到上证综合指数样本概率密度函数与基于广义双曲线分布的概率密度函 数的图形比较。由图可以看出：广义双曲线很好地拟合了样本的经验分布。 函数 gh_pdf 为广义双曲线分布的概率密度函数，其具体内容如下： function [xx,pdf]=gh_pdf(x,y) xx=sort(y); lambda=x(1);alpha=x(2);beta=x(3);mu=x(4);delta=x(5); a=(alpha^2-beta^2)^(lambda/2); b=sqrt(2*pi)*alpha.^(lambda-0.5)*delta^(lambda); c=besselk(lambda,delta*sqrt(alpha.^2-beta^2)); g=a/(b*c); a1=(delta^2+(xx-mu).^2).^(lambda/2-0.25); b1=besselk(lambda-0.5,alpha*sqrt(delta^2+(xx-mu).^2)); c1=exp(beta*(xx-mu)); pdf=g*a1.*b1.*c1;