arma模型程序

A R M A模型程序(总6页)
--本页仅作为文档封面，使用时请直接删除即可--
--内页可以根据需求调整合适字体及大小--
clear,clc
close all
set(0,'defaultaxeslinestyleorder',{'',':+','-s'}); %设置默认的线性和标识符
set(0,'defaultaxescolororder',[0,0,1])%设置默认的线条颜色
x=-pi:pi/10:pi;
y=tan(sin(x))-sin(tan(x));
set(gca,'xtick',[min(x),max(x)]) %设置x轴的范围
set(gca,'ylim',[ min(y) ,max(y)],'layer','top')%设置y轴的范围
h=plot(x,y,'markersize',9); %画图并且取得二维图形中最大值和最小值的索引号grid on
set(gca,'layer','top') %将栅格放在最上层，防止遮盖
axis equal %使每一格的刻度相同，同时也使长度一样
xlabel('-\pi\leq \theta \leq\pi','fontsize',20)
ylabel('tan(sin(\theta))-sin(tan(\theta))','fontsize',20)
title('\it{-\pi\leq \theta\leq\pi的三角函数图}','fontsize',20)
text(0,0,'\bullet\it原点','fontsize',25) %文本的精确定位
%标注曲线的最大值和最小值
x=get(h,'xdata'); %获得二维曲线的数据
y=get(h,'ydata');
imin=find(min(y)==y);
imax=find(max(y)==y);
text(x(imin),y(imin),['最小值=',num2str(y(imin))],...
'horizontalalignment','center','verticalalignment','bottom','fontsize',17) text(x(imax),y(imax),['最大值=',num2str(y(imax))],...
'horizontalalignment','center','verticalalignment','top','fontsize',17)
宋哲 2011/8/16 11:07:42
%AR,MA模型
clear,clc
close all
a=load('');
a=a([2,4],:);a=a';
a=a(:);%将数据税收展成按时间序列排列的
% plot(a,'-p')
b=diff(a);
% plot(b,'-p')
r1=autocorr(b)%计算自相关系数
r2=parcorr(b) %计算偏自相关系数
figure, subplot(121),autocorr(b)
subplot(122), parcorr(b)
cs1=ar(b,2) %AR模型，b必须是列向量
bhat=predict(cs1,[b;20]) %Discrete-time IDPOLY model:
%A(q) = 1 - q^-1 - q^-2
b15=bhat{1}(end)
bhat2=predict(cs1,[b;b15;-1])
宋哲 2011/8/16 11:27:00
%用逻辑斯特预测
clear,clc
a=load('');
a=a([2,4],:);a=a';
a=a(:);%将数据税收展成按时间序列排列的
x0=a(1);tt=[1:14]';
t0=1;
xt=@(cs,t)cs(1)/(1+(cs(1)/x0-1)*exp(-cs(2)*(t-
t0)));cs=lsqcurvefit(xt,rand(2,1),tt(2,end),a(2:end))
%以上程序为三维视图和等高线画图程序
clc, clear
a=load('');
amin=min(min(a)), amax=max(max(a))
x0=0:5:4000; y0=0:5:3000;
[xx0,yy0]=meshgrid(x0,y0);
mesh(xx0,yy0,a)
v=[-49,1,12,44,75,107,171,184:80:357,357:100:614,614];
figure, c=contour(x0,y0,a,v); clabel(c)
宋哲 2011-08-17 9:32:44
clc, clear %该程序计算最小面积
a=load('');
a=a';
x0=0:5:4000; y0=0:5:3000;
pp=csape({x0,y0},a)
x=0:4000; y=0:3000;
z=fnval(pp,{x,y});
[m,n]=size(z);
s=0;
for i=1:m-1
for j=1:n-1
p1=[x(i),y(j),z(i,j)];
p2=[x(i+1),y(j),z(i+1,j)];
p3=[x(i+1),y(j+1),z(i+1,j+1)];
p4=[x(i),y(j+1),z(i,j+1)];
p12=norm(p1-p2); p23=norm(p3-p2); p13=norm(p3-p1); p14=norm(p4-p1); p34=norm(p4-p3);
if all([z(i,j),z(i+1,j),z(i+1,j+1)]>=12)
z1=(p12+p23+p13)/2;
s1=sqrt(z1*(z1-p12)*(z1-p23)*(z1-p13));
s=s+s1;
end
if all([z(i,j),z(i+1,j+1),z(i,j+1)]>=12)
z2=(p13+p14+p34)/2;
s2=sqrt(z2*(z2-p13)*(z2-p14)*(z2-p34));
s=s+s2;
end
end
end
s
宋哲 2011-08-17 9:32:57
%图3的Matlab程序如下
clc, clear
load pdata
a=load(' ');
x0=0:5:4000; y0=0:5:3000;
v=[-49,1,12,44,75,107,171,184:80:357,357:100:614,614];
c=contour(x0,y0,a,v); clabel(c)
hold on
plot(P1(1),P1(2),'S')
plot(P2(1),P2(2),'D'), plot(P3(1),P3(2),'D'), plot(P4(1),P4(2),'D') b=a';
for i=1:801
for j=1:601
txy=[(i-1)*5;(j-1)*5;b(i,j)]; td=norm(txy-Pz1);
if td>=995 & td<=1000 plot((i-1)*5,(j-1)*5,'o') end
end
end
for i=1:801
for j=1:601
txy=[(i-1)*5;(j-1)*5;b(i,j)]; td=norm(txy-Pz2);
if td>=495 & td<=500 plot((i-1)*5,(j-1)*5,'o') end
end
end
for i=1:801
for j=1:601
txy=[(i-1)*5;(j-1)*5;b(i,j)]; td=norm(txy-Pz3);
if td>=495 & td<=500
plot((i-1)*5,(j-1)*5,'o')
end
end
end
for i=1:801
for j=1:601
txy=[(i-1)*5;(j-1)*5;b(i,j)];
td=norm(txy-Pz4);
if td>=495 & td<=500
plot((i-1)*5,(j-1)*5,'o')
end
end
end
宋哲 2011-08-17 16:11:05
clc, clear %该程序计算最小面积a=load('');
a=a';
x0=0:5:4000; y0=0:5:3000;
pp=csape({x0,y0},a)
x=0:4000; y=0:3000;
z=fnval(pp,{x,y});
[m,n]=size(z)
s=0;
for i=1:m-1
for j=1:n-1
p1=[x(i),y(j),z(i,j)];
p2=[x(i+1),y(j),z(i+1,j)];
p3=[x(i+1),y(j+1),z(i+1,j+1)];
p4=[x(i),y(j+1),z(i,j+1)];
p12=norm(p1-p2); p23=norm(p3-p2); p13=norm(p3-p1); p14=norm(p4-p1); p34=norm(p4-p3);
if any([z(i,j),z(i+1,j),z(i+1,j+1)]>=12)
z1=(p12+p23+p13)/2;
s1=sqrt(z1*(z1-p12)*(z1-p23)*(z1-p13));
s=s+s1;
end
if any([z(i,j),z(i+1,j+1),z(i,j+1)]>=12)
z2=(p13+p14+p34)/2;
s2=sqrt(z2*(z2-p13)*(z2-p14)*(z2-p34));
s=s+s2;
end
end
end
s
宋哲 2011-08-17 16:41:37
%求4个特殊点的Matlab程序
clc, clear
a=load('');
a=a';
x0=0:5:4000; y0=0:5:3000;
i=[2500:5:3700]/5; j=[800:5:2100]/5; %求358高地的大致范围的地址z0=a(i,j); %求358高地的大致高程
zz1=max(max(z0)) %求358高地的最高点高程
[i1,j1]=find(a==zz1); %求358高地在全部数据矩阵a中的地址
P1=[(i1-1)*5; (j1-1)*5] %求358高地的x坐标和y坐标
zz2=max(max(a))
[i2,j2]=find(a==zz2);
P2=[(i2-1)*5; (j2-1)*5] %求615高地的最高点x坐标和y坐标
a0=a([800:5:1100]/5,[600:5:800]/5);
zz3=max(max(a0)) %求185高地的最高点高程
[i3,j3]=find(a==zz3);
P3=[(i3-1)*5; (j3-1)*5] %求185高地的最高点x坐标和y坐标
z4=a([2100:5:2650]/5,[170:5:580]/5);
zz4=max(max(z4))
[i4,j4]=find(a==zz4);
P4=[(i4-1)*5; (j4-1)*5]
Pz1=[P1;zz1+10]; Pz2=[P2; zz2+3]; Pz3=[P3; zz3+3]; Pz4=[P4; zz4+3]; save pdata P1 P2 P3 P4 Pz1 Pz2 Pz3 Pz4
宋哲 2011-08-17 16:41:48
%图3的Matlab程序如下
clc, clear
load pdata
a=load('');
x0=0:5:4000; y0=0:5:3000;
v=[-49,1,12,44,75,107,171,184:80:357,357:100:614,614];
c=contour(x0,y0,a,v); clabel(c)
hold on
plot(P1(1),P1(2),'S')
plot(P2(1),P2(2),'D'), plot(P3(1),P3(2),'D'), plot(P4(1),P4(2),'D')
b=a';
for i=1:801
for j=1:601
txy=[(i-1)*5;(j-1)*5;b(i,j)];
td=norm(txy-Pz1);
if td>=995 & td<=1000
plot((i-1)*5,(j-1)*5,'o')
end
end
end
for i=1:801
for j=1:601
txy=[(i-1)*5;(j-1)*5;b(i,j)]; td=norm(txy-Pz2);
if td>=495 & td<=500 plot((i-1)*5,(j-1)*5,'o') end
end
end
for i=1:801
for j=1:601
txy=[(i-1)*5;(j-1)*5;b(i,j)]; td=norm(txy-Pz3);
if td>=495 & td<=500 plot((i-1)*5,(j-1)*5,'o') end
end
end
for i=1:801
for j=1:601
txy=[(i-1)*5;(j-1)*5;b(i,j)]; td=norm(txy-Pz4);
if td>=495 & td<=500 plot((i-1)*5,(j-1)*5,'o') end
end
end。