%===============================================================
clear;
syms rou fai2 k1 k2 k3 n rorn ii
clc;
n=input('观测时刻数 n=');
disp('= = = = = = = = = = = = dealing = = = = = = = = = = = = = = ');
disp('Just wait for a few minutes............');
k1=sym('(1-rou*rou)*(1-fai2*fai2)*Xmn(rorn,1)*Xmn(rorn,1)');
k2=sym('(1-fai2^2)*(Xmn(rorn,2)-rou*Xmn(rorn,1))^2');
% digits(25); %设定计算精度
k3=sym('(Xmn(rorn,ii)-rou*(1-fai2)*Xmn(rorn,ii-1)-fai2*Xmn(rorn,ii-2))^2');
k3=symsum(k3,ii,3,n);
sigema=1/n*(k1+k2+k3); %1.5
sigema_rou=1/n*(diff(k1,rou)+diff(k2,rou)+diff(k3,rou));
sigema_fai2=1/n*(diff(k1,fai2)+diff(k2,fai2)+diff(k3,fai2));
f1=(n/(2*sigema))*sigema_rou-rou/(1-rou*rou); %1.7
f2=(n/(2*sigema))*sigema_fai2-2*fai2/(1-fai2*fai2);
disp('= = = = = = = = = = = = funcation F = = = = = = = = = = = = = = ');
f1=subs(f1,'rou_fai2(1)','rou');
f2=subs(f2,'rou_fai2(1)','rou');
f1=subs(f1,'rou_fai2(2)','fai2');
f2=subs(f2,'rou_fai2(2)','fai2');
F=[f1;f2];
disp('F=');
disp(F);
disp('= = = = = = = = = = = = simpling...... = = = = = = = = = = = = = = ');
f1=simple(f1);
f2=simple(f2);
F=[f1;f2];
disp('F=');
disp(F);
%===============================================================
second_step.m
%===============================================================
%极大近似估计函数
%个参数在程序中的替代表示：δ=sigema, ξ=ipsn,φ=fai,ρ=rou
%note:sigema在程序中实际表示的是δ的平方
% saved_max_approx.m的结果填入F中。
function F=second_step(rou_fai2) %max approximate
global rorn Xmn
%*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*
%copy the value of F in the first_step to fill in F=[] as following:
%for instance,n=10:
% F=[10/(1/5*(1-rou_fai2(1)^2)*(1-rou_fai2(2)^2)*Xmn(rorn,1)^2+1/5*(1-rou_fai2(2)^2)*(Xmn(rorn,2)-rou_fai2(1)*Xmn(rorn,1))^2+1/5*(Xmn(rorn,3)-rou_fai2(1)*(1-rou_fai2(2))*Xmn(rorn,2)-rou_fai2(2)*Xmn(rorn,1))^2+1/5*(Xmn(rorn,4)-rou_fai2(1)*(1-rou_fai2(2))*Xmn(rorn,3)-rou_fai2(2)*Xmn(rorn,2))^2+1/5*(Xmn(rorn,5)-rou_fai2(1)*(1-rou_fai2(2))*Xmn(rorn,4)-rou_fai2(2)*Xmn(rorn,3))^2+1/5*(Xmn(rorn,6)-rou_fai2(1)*(1-rou_fai2(2))*Xmn(rorn,5)-rou_fai2(2)*Xmn(rorn,4))^2+1/5*(Xmn(rorn,7)-rou_fai2(1)*(1-rou_fai2(2))*Xmn(rorn,6)-rou_fai2(2)*Xmn(rorn,5))^2+1/5*(Xmn(rorn,8)-rou_fai2(1)*(1-rou_fai2(2))*Xmn(rorn,7)-rou_fai2(2)*Xmn(rorn,6))^2+1/5*(Xmn(rorn,9)-rou_fai2(1)*(1-rou_fai2(2))*Xmn(rorn,8)-rou_fai2(2)*Xmn(rorn,7))^2+1/5*(Xmn(rorn,10)-rou_fai2(1)*(1-rou_fai2(2))*Xmn(rorn,9)-rou_fai2(2)*Xmn(rorn,8))^2)*(-1/5*rou_fai2(1)*(1-rou_fai2(2)^2)*Xmn(rorn,1)^2-1/5*(1-rou_fai2(2)^2)*(Xmn(rorn,2)-rou_fai2(1)*Xmn(rorn,1))*Xmn(rorn,1)-1/5*(Xmn(rorn,3)-rou_fai2(1)*(1-rou_fai2(2))*Xmn(rorn,2)-rou_fai2(2)*Xmn(rorn,1))*(1-rou_fai2(2))*Xmn(rorn,2)-1/5*(Xmn(rorn,4)-rou_fai2(1)*(1-rou_fai2(2))*Xmn(rorn,3)-rou_fai2(2)*Xmn(rorn,2))*(1-rou_fai2(2))*Xmn(rorn,3)-1/5*(Xmn(rorn,5)-rou_fai2(1)*(1-rou_fai2(2))*Xmn(rorn,4)-rou_fai2(2)*Xmn(rorn,3))*(1-rou_fai2(2))*Xmn(rorn,4)-1/5*(Xmn(rorn,6)-rou_fai2(1)*(1-rou_fai2(2))*Xmn(rorn,5)-rou_fai2(2)*Xmn(rorn,4))*(1-rou_fai2(2))*Xmn(rorn,5)-1/5*(Xmn(rorn,7)-rou_fai2(1)*(1-rou_fai2(2))*Xmn(rorn,6)-rou_fai2(2)*Xmn(rorn,5))*(1-rou_fai2(2))*Xmn(rorn,6)-1/5*(Xmn(rorn,8)-rou_fai2(1)*(1-rou_fai2(2))*Xmn(rorn,7)-rou_fai2(2)*Xmn(rorn,6))*(1-rou_fai2(2))*Xmn(rorn,7)-1/5*(Xmn(rorn,9)-rou_fai2(1)*(1-rou_fai2(2))*Xmn(rorn,8)-rou_fai2(2)*Xmn(rorn,7))*(1-rou_fai2(2))*Xmn(rorn,8)-1/5*(Xmn(rorn,10)-rou_fai2(1)*(1-rou_fai2(2))*Xmn(rorn,9)-rou_fai2(2)*Xmn(rorn,8))*(1-rou_fai2(2))*Xmn(rorn,9))-rou_fai2(1)/(1-rou_fai2(1)^2)
% 10/(1/5*(1-rou_fai2(1)^2)*(1-rou_fai2(2)^2)*Xmn(rorn,1)^2+1/5*(1-rou_fai2(2)^2)*(Xmn(rorn,2)-rou_fai2(1)*Xmn(rorn,1))^2+1/5*(Xmn(rorn,3)-rou_fai2(1)*(1-rou_fai2(2))*Xmn(rorn,2)-rou_fai2(2)*Xmn(rorn,1))^2+1/5*(Xmn(rorn,4)-rou_fai2(1)*(1-rou_fai2(2))*Xmn(rorn,3)-rou_fai2(2)*Xmn(rorn,2))^2+1/5*(Xmn(rorn,5)-rou_fai2(1)*(1-rou_fai2(2))*Xmn(rorn,4)-rou_fai2(2)*Xmn(rorn,3))^2+1/5*(Xmn(rorn,6)-rou_fai2(1)*(1-rou_fai2(2))*Xmn(rorn,5)-rou_fai2(2)*Xmn(rorn,4))^2+1/5*(Xmn(rorn,7)-rou_fai2(1)*(1-rou_fai2(2))*Xmn(rorn,6)-rou_fai2(2)*Xmn(rorn,5))^2+1/5*(Xmn(rorn,8)-rou_fai2(1)*(1-rou_fai2(2))*Xmn(rorn,7)-rou_fai2(2)*Xmn(rorn,6))^2+1/5*(Xmn(rorn,9)-rou_fai2(1)*(1-rou_fai2(2))*Xmn(rorn,8)-rou_fai2(2)*Xmn(rorn,7))^2+1/5*(Xmn(rorn,10)-rou_fai2(1)*(1-rou_fai2(2))*Xmn(rorn,9)-rou_fai2(2)*Xmn(rorn,8))^2)*(-1/5*(1-rou_fai2(1)^2)*rou_fai2(2)*Xmn(rorn,1)^2-1/5*rou_fai2(2)*(Xmn(rorn,2)-rou_fai2(1)*Xmn(rorn,1))^2+1/5*(Xmn(rorn,3)-rou_fai2(1)*(1-rou_fai2(2))*Xmn(rorn,2)-rou_fai2(2)*Xmn(rorn,1))*(rou_fai2(1)*Xmn(rorn,2)-Xmn(rorn,1))+1/5*(Xmn(rorn,4)-rou_fai2(1)*(1-rou_fai2(2))*Xmn(rorn,3)-rou_fai2(2)*Xmn(rorn,2))*(rou_fai2(1)*Xmn(rorn,3)-Xmn(rorn,2))+1/5*(Xmn(rorn,5)-rou_fai2(1)*(1-rou_fai2(2))*Xmn(rorn,4)-rou_fai2(2)*Xmn(rorn,3))*(rou_fai2(1)*Xmn(rorn,4)-Xmn(rorn,3))+1/5*(Xmn(rorn,6)-rou_fai2(1)*(1-rou_fai2(2))*Xmn(rorn,5)-rou_fai2(2)*Xmn(rorn,4))*(rou_fai2(1)*Xmn(rorn,5)-Xmn(rorn,4))+1/5*(Xmn(rorn,7)-rou_fai2(1)*(1-rou_fai2(2))*Xmn(rorn,6)-rou_fai2(2)*Xmn(rorn,5))*(rou_fai2(1)*Xmn(rorn,6)-Xmn(rorn,5))+1/5*(Xmn(rorn,8)-rou_fai2(1)*(1-rou_fai2(2))*Xmn(rorn,7)-rou_fai2(2)*Xmn(rorn,6))*(rou_fai2(1)*Xmn(rorn,7)-Xmn(rorn,6))+1/5*(Xmn(rorn,9)-rou_fai2(1)*(1-rou_fai2(2))*Xmn(rorn,8)-rou_fai2(2)*Xmn(rorn,7))*(rou_fai2(1)*Xmn(rorn,8)-Xmn(rorn,7))+1/5*(Xmn(rorn,10)-rou_fai2(1)*(1-rou_fai2(2))*Xmn(rorn,9)-rou_fai2(2)*Xmn(rorn,8))*(rou_fai2(1)*Xmn(rorn,9)-Xmn(rorn,8)))-2*rou_fai2(2)/(1-rou_fai2(2)^2)];
F=[10/(1/5*(1-(rou_fai2(1))^2)*(1-(rou_fai2(2))^2)*Xmn(rorn,1)^2+1/5*(1-(rou_fai2(2))^2)*(Xmn(rorn,2)-(rou_fai2(1))*Xmn(rorn,1))^2+1/5*(Xmn(rorn,3)-(rou_fai2(1))*(1-(rou_fai2(2)))*Xmn(rorn,2)-(rou_fai2(2))*Xmn(rorn,1))^2+1/5*(Xmn(rorn,4)-(rou_fai2(1))*(1-(rou_fai2(2)))*Xmn(rorn,3)-(rou_fai2(2))*Xmn(rorn,2))^2+1/5*(Xmn(rorn,5)-(rou_fai2(1))*(1-(rou_fai2(2)))*Xmn(rorn,4)-(rou_fai2(2))*Xmn(rorn,3))^2+1/5*(Xmn(rorn,6)-(rou_fai2(1))*(1-(rou_fai2(2)))*Xmn(rorn,5)-(rou_fai2(2))*Xmn(rorn,4))^2+1/5*(Xmn(rorn,7)-(rou_fai2(1))*(1-(rou_fai2(2)))*Xmn(rorn,6)-(rou_fai2(2))*Xmn(rorn,5))^2+1/5*(Xmn(rorn,8)-(rou_fai2(1))*(1-(rou_fai2(2)))*Xmn(rorn,7)-(rou_fai2(2))*Xmn(rorn,6))^2+1/5*(Xmn(rorn,9)-(rou_fai2(1))*(1-(rou_fai2(2)))*Xmn(rorn,8)-(rou_fai2(2))*Xmn(rorn,7))^2+1/5*(Xmn(rorn,10)-(rou_fai2(1))*(1-(rou_fai2(2)))*Xmn(rorn,9)-(rou_fai2(2))*Xmn(rorn,8))^2)*(-1/5*(rou_fai2(1))*(1-(rou_fai2(2))^2)*Xmn(rorn,1)^2-1/5*(1-(rou_fai2(2))^2)*(Xmn(rorn,2)-(rou_fai2(1))*Xmn(rorn,1))*Xmn(rorn,1)-1/5*(Xmn(rorn,3)-(rou_fai2(1))*(1-(rou_fai2(2)))*Xmn(rorn,2)-(rou_fai2(2))*Xmn(rorn,1))*(1-(rou_fai2(2)))*Xmn(rorn,2)-1/5*(Xmn(rorn,4)-(rou_fai2(1))*(1-(rou_fai2(2)))*Xmn(rorn,3)-(rou_fai2(2))*Xmn(rorn,2))*(1-(rou_fai2(2)))*Xmn(rorn,3)-1/5*(Xmn(rorn,5)-(rou_fai2(1))*(1-(rou_fai2(2)))*Xmn(rorn,4)-(rou_fai2(2))*Xmn(rorn,3))*(1-(rou_fai2(2)))*Xmn(rorn,4)-1/5*(Xmn(rorn,6)-(rou_fai2(1))*(1-(rou_fai2(2)))*Xmn(rorn,5)-(rou_fai2(2))*Xmn(rorn,4))*(1-(rou_fai2(2)))*Xmn(rorn,5)-1/5*(Xmn(rorn,7)-(rou_fai2(1))*(1-(rou_fai2(2)))*Xmn(rorn,6)-(rou_fai2(2))*Xmn(rorn,5))*(1-(rou_fai2(2)))*Xmn(rorn,6)-1/5*(Xmn(rorn,8)-(rou_fai2(1))*(1-(rou_fai2(2)))*Xmn(rorn,7)-(rou_fai2(2))*Xmn(rorn,6))*(1-(rou_fai2(2)))*Xmn(rorn,7)-1/5*(Xmn(rorn,9)-(rou_fai2(1))*(1-(rou_fai2(2)))*Xmn(rorn,8)-(rou_fai2(2))*Xmn(rorn,7))*(1-(rou_fai2(2)))*Xmn(rorn,8)-1/5*(Xmn(rorn,10)-(rou_fai2(1))*(1-(rou_fai2(2)))*Xmn(rorn,9)-(rou_fai2(2))*Xmn(rorn,8))*(1-(rou_fai2(2)))*Xmn(rorn,9))-(rou_fai2(1))/(1-(rou_fai2(1))^2)
10/(1/5*(1-(rou_fai2(1))^2)*(1-(rou_fai2(2))^2)*Xmn(rorn,1)^2+1/5*(1-(rou_fai2(2))^2)*(Xmn(rorn,2)-(rou_fai2(1))*Xmn(rorn,1))^2+1/5*(Xmn(rorn,3)-(rou_fai2(1))*(1-(rou_fai2(2)))*Xmn(rorn,2)-(rou_fai2(2))*Xmn(rorn,1))^2+1/5*(Xmn(rorn,4)-(rou_fai2(1))*(1-(rou_fai2(2)))*Xmn(rorn,3)-(rou_fai2(2))*Xmn(rorn,2))^2+1/5*(Xmn(rorn,5)-(rou_fai2(1))*(1-(rou_fai2(2)))*Xmn(rorn,4)-(rou_fai2(2))*Xmn(rorn,3))^2+1/5*(Xmn(rorn,6)-(rou_fai2(1))*(1-(rou_fai2(2)))*Xmn(rorn,5)-(rou_fai2(2))*Xmn(rorn,4))^2+1/5*(Xmn(rorn,7)-(rou_fai2(1))*(1-(rou_fai2(2)))*Xmn(rorn,6)-(rou_fai2(2))*Xmn(rorn,5))^2+1/5*(Xmn(rorn,8)-(rou_fai2(1))*(1-(rou_fai2(2)))*Xmn(rorn,7)-(rou_fai2(2))*Xmn(rorn,6))^2+1/5*(Xmn(rorn,9)-(rou_fai2(1))*(1-(rou_fai2(2)))*Xmn(rorn,8)-(rou_fai2(2))*Xmn(rorn,7))^2+1/5*(Xmn(rorn,10)-(rou_fai2(1))*(1-(rou_fai2(2)))*Xmn(rorn,9)-(rou_fai2(2))*Xmn(rorn,8))^2)*(-1/5*(1-(rou_fai2(1))^2)*(rou_fai2(2))*Xmn(rorn,1)^2-1/5*(rou_fai2(2))*(Xmn(rorn,2)-(rou_fai2(1))*Xmn(rorn,1))^2+1/5*(Xmn(rorn,3)-(rou_fai2(1))*(1-(rou_fai2(2)))*Xmn(rorn,2)-(rou_fai2(2))*Xmn(rorn,1))*((rou_fai2(1))*Xmn(rorn,2)-Xmn(rorn,1))+1/5*(Xmn(rorn,4)-(rou_fai2(1))*(1-(rou_fai2(2)))*Xmn(rorn,3)-(rou_fai2(2))*Xmn(rorn,2))*((rou_fai2(1))*Xmn(rorn,3)-Xmn(rorn,2))+1/5*(Xmn(rorn,5)-(rou_fai2(1))*(1-(rou_fai2(2)))*Xmn(rorn,4)-(rou_fai2(2))*Xmn(rorn,3))*((rou_fai2(1))*Xmn(rorn,4)-Xmn(rorn,3))+1/5*(Xmn(rorn,6)-(rou_fai2(1))*(1-(rou_fai2(2)))*Xmn(rorn,5)-(rou_fai2(2))*Xmn(rorn,4))*((rou_fai2(1))*Xmn(rorn,5)-Xmn(rorn,4))+1/5*(Xmn(rorn,7)-(rou_fai2(1))*(1-(rou_fai2(2)))*Xmn(rorn,6)-(rou_fai2(2))*Xmn(rorn,5))*((rou_fai2(1))*Xmn(rorn,6)-Xmn(rorn,5))+1/5*(Xmn(rorn,8)-(rou_fai2(1))*(1-(rou_fai2(2)))*Xmn(rorn,7)-(rou_fai2(2))*Xmn(rorn,6))*((rou_fai2(1))*Xmn(rorn,7)-Xmn(rorn,6))+1/5*(Xmn(rorn,9)-(rou_fai2(1))*(1-(rou_fai2(2)))*Xmn(rorn,8)-(rou_fai2(2))*Xmn(rorn,7))*((rou_fai2(1))*Xmn(rorn,8)-Xmn(rorn,7))+1/5*(Xmn(rorn,10)-(rou_fai2(1))*(1-(rou_fai2(2)))*Xmn(rorn,9)-(rou_fai2(2))*Xmn(rorn,8))*((rou_fai2(1))*Xmn(rorn,9)-Xmn(rorn,8)))-2*(rou_fai2(2))/(1-(rou_fai2(2))^2)];
%It's a function,so needn't execute
%===============================================================
third_step.m
%===============================================================
%矩阵运算之解方程组
%输入参数：n为整型数，φ1、φ2，ε1～εn为双精度浮点数，须带小数点的
%输出数据：X1～Xn
%================================================================
clc;
disp('! 注意：是否将first_step中产生的F的两列表达式复制到程序second_step中且保存！请选择y/n:');
userchoice = input('@=@=@=@=@=@=@=@=@=@=@=@=选择：','s');
if userchoice ~='y';
error('请重新运行程序first_step，并将结果填入程序second_step中F=[]处！');
end
clc;
disp(' 输入预设数据：');%disp(['',n,'']);
% n=input('输入参数n='); %n is filled in the first_step
m=input('输入参数m=');
p=input('输入参数φ1=');
q=input('输入参数φ2=');
si=input('输入参数si=');%si代表normrnd函数的参数sigema
%用normrnd（0，sigma，行，列）产生正态随机数
%矩阵下标从1开始
%================================================================
E=normrnd(0,si,m,n); %按序产生随机数ε1～εn,m*n矩阵
disp('E(m,n)=');disp(E);
X=zeros(m,n);
Xn=zeros(2,n);
for k=1:m
for i0=1:n
if i0==1;
Xn(2,1)=E(k,1);
X1=Xn(2,1);
end
if i0==2;
Xn(2,2)=p*Xn(2,1)+E(k,2);
X2=Xn(2,2);
end
if i0>2;
Xn(2,i0)=p*Xn(2,i0-1)+q*Xn(2,i0-2)+E(k,i0);
end
if i0==n;
for i1=1:n;
X(k,i1)=Xn(2,i1);
end
end
% i0=i0+1;?
end
end
disp('X=');
disp(X);
disp('X dealing finished!');
%================================================================
%warning('message');%不退出程序
%error('message'); %退出程序
%break; %跳出for或while循环
%================================================================
%下面用fsolve函数来解方程组1.8
rou=zeros(1,m);
fai2=zeros(1,m);
global rorn Xmn
Xmn=X;
options=optimset('Display','iter'); % Option to display output
rou_fai2_0 = [0.5; 0.5]; %起始数据，可作为调整参数
for rorn0=1:m;
rorn=rorn0;
disp('*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*');
disp('当前程序执行到的X矩阵的行数为：');disp(rorn);
[rou_fai2,fval] = fsolve(@second_step,rou_fai2_0,options) ;
disp('rou=');disp(rou_fai2(1));rou(rorn)=rou_fai2(1);
disp('fai2=');disp(rou_fai2(2));fai2(rorn)=rou_fai2(2);
disp('fval=');disp(fval);
end
%满足极大似然估计的（rou,fai2)的判定，以及后续处理：
sigema_2=zeros(1,m);
fai_1=zeros(1,m);
for rorm=1:m;
disp('=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=');
disp('数据在执行到行数为');disp(rorm);disp(' 时，有处理结果为：')
if abs(rou(rorm))<1;
if abs(fai2(rorm))<1;
sigema2=sigema; %1.9
sigema2=subs(sigema2,'rou(rorm)','rou');
sigema2=subs(sigema2,'fai2(rorm)','fai2');
sigema2=simple(sigema2);
sigema2=eval(sigema2);
fai1=rou(rorm)*(1-fai2(rorm)); %1.10
sigema_2(rorm)=sigema2;
fai_1(rorm)=fai1;
disp('有满足条件的结果！！ 结果为：')
disp('σ^2=');disp(sigema2);
disp('φ1=');disp(fai1);
else
disp('|φ2|>=1,不满足集合条件！');
end
else
disp('|ρ|>=1,不满足集合条件！');
end
end
disp('=*=*=*=*=*=*=*=*=*=*=*The result!*=*=*=*=*=*=*=*=*=*=');
disp('rou=');disp(rou);
disp('φ2=');disp(fai2);
disp('σ^2的极大似然估计为');
disp(sigema_2);
disp('φ1的极大似然估计为');
disp(fai_1);
%将结果输出到excel：
xlswrite('X', X);
xlswrite('rou.xls', rou);
xlswrite('fai2.xls', fai2);
xlswrite('sigema_2.xls', sigema_2);
xlswrite('fai1', fai_1);
% var %求方差
% std %求标准差
% % 绘制近似数图像
% title('φ1、φ2的散点图');grid;
% subplot(2,1,1);
% len_fai=length(fai_1);
% n0=1:len_fai;
% scatter(n0,fai_1);
% hold on;
% n1=fix(len_fai/2);
% stem(n1,p);
%
% subplot(2,1,2);
% scatter(n0,fai2);
% hold on;
% stem(n1,q);
%===============================================================
