VASICEK模型GMM估计时出现问题，真心求助！-经管之家官网！

我参考了网上曹老师的程序后，自己在细节的地方修改了一下，但是不知道哪里出现了错误，一直提示错误，希望懂的同学和老师可以帮帮忙，先谢谢大家了！
function [paraest,t_sta,V,it,Chi_sta,Pvalue]=gmmestimation
(moment,para0,Y,X,Z,number,K)
nlag=round(size(Y,1)^(1/3));
W(:,:,1)=eye(K);
[para(:,1),fv(:,1)]=fminsearch(moment,para(:,1),[],1,Y,X,Z,W(:,:,1));
% update weight matrix and find optimal solution
for i=2:number
mom=feval_r(moment,para(:,i-1),2,Y,X,Z,W(:,:,i-1));
W(:,:,i)=gmmweightmatrix(mom,nlag);
[para(:,i),fv(:,i)]=fminsearch(moment,para(:,i-1),[],1,Y,X,Z,W(:,:,i));
if abs(fv(:,i)-fv(:,i-1))/abs(fv(:,i-1))<1e-4|fv(i)<=1e-10
break
end
end
it=i;
if it==number
error('number of iteration exceeds defined maximum number')
else
paraest=para(:,it);% optimal parameter
f0=feval_r(moment,paraest,3,Y,X,Z,W(:,:,it));% optimal function value
% find covariance matrix of estimated parameter vector, using
numerical method
for j=1:length(para0)
a=zeros(length(para0),1);
eps=max(paraest(j)*1e-6,1e-5);
a(j)=eps;
M(:,j)=(feval_r(moment,paraest+a,3,Y,X,Z,W(:,:,it))-f0)/eps;
end
end
V=pinv(M'*W(:,:,it)*M)/size(Y,1);
stderror=sqrt(diag(V));
t_sta=paraest./(stderror);
Chi_sta=size(Y,1)*fv(it);
Pvalue=1-chi2cdf(Chi_sta,K-length(para0));
%find the weight matrix using Newey and West method
function W=gmmweightmatrix(mom,nlag)
q=size(mom,2);T=size(mom,1);
a2=zeros(q,q);a3=zeros(q,q);
for j=1:nlag
a1=zeros(q,q);
for i=1:(T-j)
a1=mom(i+j,:)'*mom(i,:)+a1;
end
S(:,:,j)=1/T*a1;
a2=(1-j/(nlag+1))*S(:,:,j)+a2;
a3=(1-j/(nlag+1))*S(:,:,j)'+a3;
end
b1=zeros(q,q);
for i=1:T
b1=mom(i,:)'*mom(i,:)+b1;
end
if nlag==0
newS=b1*1/T;
else
newS=a2+a3+b1*1/T;
end
W=pinv(newS);
%%%%%%%%%%%%%%%%%%%%%
function f=cklstest01(para,num,Y,X,Z,W)
[T,q]=size(Y);
alpha=para(1);beta=para(2);sigsq=para(3);
eta=[Y-(alpha+beta*X)/12,(Y-(alpha+beta*X)/12).^2-sigsq.^2/12];
for i=1:T
m_t(i,:)=kron(eta(i,:),Z(i,:));
end
m=mean(m_t)';
obj=m'*W*m;
if num==1
f=obj;
elseif num==2
f=m_t;
elseif num==3
f=m;
end
%%%%%%%%%%%%%
clear
data=xlsread('gmm.xlsx');
Y=data(:,1);
X=data(:,2);
para0=[0.05;-0.6;0.5];number=20;moment='cklstest01';K=3;
Z=[ones(length(Y),1),X];
[beta,stderror,V,it,Chi_sta,Pvalue]=gmmestimation
(moment,para0,Y,X,Z,number,K);