我做的过程是：
将y.xls导入matlab
将skewtdis_LL.m, my_mle.m加入MATLAB.然后在命令窗口输入
[para,standard_deviation]=my_mle('skewtdis_LL',[5;0.02],y)
结果显示：
??? Error using ==> skewtdis_LL
Too many input arguments.
Error in ==> fminsearch at 195
fv(:,1) = funfcn(x,varargin{:});
Error in ==> my_mle at 3
[para,fv]=fminsearch(fun,para0,[],2,varargin{:});


程序如下：
function LL = skewtdis_LL(theta,x)
% PURPOSE: returns the log-likelihood at x of Hansen's (1994) 'skewed t' distribution
%---------------------------------------------------
% USAGE: LL = skewtdis_LL(x,nu,lambda)
% where: x  = a vector of data
%        theta = [nu;lambda]
%     nu = degrees of freedom parameter
%       lambda = skewness parameter
%---------------------------------------------------
% RETURNS:
%        a matrix of log-likelihood at each element of x
% --------------------------------------------------
% SEE ALSO: tdis_cdf, tdis_rnd, tdis_inv, tdis_prb, skewtdis_pdf
%---------------------------------------------------
%
%  Andrew Patton
%
%  25 June, 2001
% This code was used in:
%
% Patton, Andrew J., 2002, "On the Out-of-Sample
% Importance of Skewness and Asymmetric Dependence for
% Asset Allocation", working paper, Department of Economics,
% University of California, San Diego.
nu = theta(1);
lambda = theta(2);
[T,k] = size(x);
nu = nu(1)*ones(T,1);     % can make this time-varying, but needs to be >2
lambda = lambda(1)*ones(T,1); % can make this time-varying, but needs to be in (-1,1)
   
c = gamma((nu+1)/2)./(sqrt(pi*(nu-2)).*gamma(nu/2));
a = 4*lambda.*c.*((nu-2)./(nu-1));
b = sqrt(1 + 3*lambda.^2 - a.^2);
logc = gammaln((nu+1)/2) - gammaln(nu/2) - 0.5*log(pi*(nu-2));
logb = 0.5*log(1 + 3*lambda.^2 - a.^2);
find1 = (x<(-a./b));
find2 = (x>=(-a./b));
LL1   = logb + logc - (nu+1)/2.*log(1+1./(nu-2).*((b.*x+a)./(1-lambda)).^2);
LL2   = logb + logc - (nu+1)/2.*log(1+1./(nu-2).*((b.*x+a)./(1+lambda)).^2);
LL    = sum(LL1(find1)) + sum(LL2(find2));
LL    = -LL;
function [para,standard_deviation,fv]=my_mle(fun,para0,varargin)
para0=para0(:);
[para,fv]=fminsearch(fun,para0,[],2,varargin{:});
fv=-fv;
d=numericalfirstderivative(fun,para,1,varargin{:});
standard_deviation=sqrt(diag(pinv(d'*d)));
function f=numericalfirstderivative(fun,parameter,varargin)
n=length(parameter);
for i=1:n
a=zeros(n,1);
a(i)=min(parameter(i)*1e-6,1e-5);
y1(:,i)=feval_r(fun,parameter+a,varargin{:});
y2(:,i)=feval_r(fun,parameter-a,varargin{:});
f(:,i)=(y1(:,i)-y2(:,i))/2/a(i);
end