function f=gm4(x,K)


n=length(x);

z(n-3,n)=0;
zz(n-2,n)=0;%
AB=zeros(2,n-3);
for j=1:n-3
     
z(j,j:n)=x(j:n);

x1=zeros(1,n-j+1);
x0=zeros(1,n-j+1);
Y=zeros(1,n-j);
B=zeros(n-j,2);%




    
x1(1)=z(j,j); 
for i=j+1:n
x1(i-j+1)=x1(i-j)+z(j,i);
Y(i-j)=x(i);
end





for i=j:n-1
B(i-j+1,1)=-0.5*(x1(i-j+1)+x1(i-j+2));
B(i-j+1,2)=1;
end






ab=inv(B'*B)*B'*Y';
a=ab(1);
b=ab(2);%
AB(1,j)=a;
AB(2,j)=b;



for i=j-1:n-1
x1(i-j+2)=(z(j,j)-b/a)*exp(-a*(i-j+1))+b/a;%
end


x0(1)=x(j);
for i=j+1:n%
x0(i-j+1)=x1(i-j+1)-x1(i-j);
end



zz(j,j:end)=x0;
end




for j=1:n
    if j<=n-3
 zz(n-2,j)=sum(zz(:,j))/j;
else
    zz(n-2,j)=sum(zz(:,j))/(n-3);
end
 
end







for j=1:n
e(j)=zz(n-2,j)-x(j);
k(j)=abs(x(j)-zz(n-2,j))/x(j);
end


'絕對誤差平方和';
s=e*e';%絕對誤差平方和
'平均相對誤差';
ss=sum(k(2:n))/(n-1);%平均相對誤差



%模擬方程
'x1(k)=(x(1)-b/a)*exp(-a*(k-1))+b/a';

'x0(k)=x1(k)-x1(k-1)';


zzz(n-3,3)=0;%第K年預測值
%預測
if K>n
 for j=1:n-3
    zzz(j,1)=(zz(j,j)-AB(2,j)/AB(1,j))*exp(-AB(1,j)*(K-j))+AB(2,j)/AB(1,j);%此處的K+1還是K一直不太清楚
    
        zzz(j,2)=(zz(j,j)-AB(2,j)/AB(1,j))*exp(-AB(1,j)*(K-j-1))+AB(2,j)/AB(1,j);
        
               zzz(j,3)= zzz(j,1)-zzz(j,2);
               zzz;
           end
           '第K年的預測值為'
           mean(zzz(:,3))
           else
            '第K年的預測值為'    
           zz(n-2,K)
       end