% 本程序的目的是模拟一个ARMA模型，然后进行时频归并。考察归并前后模型的变化。
% 这个ARMA模型的一般形式用黑盒子模型表示为A(q)y(t)=C(q)e(t)。q是滞后算子。
% 或者是：(1+a1*q^(-1)+a2*q^(-2)+a3^(-3)+a4*q^(-4))y(t)=(1+c1*q^(-1)+c2*q^(-2)+c3^(-3)+c4*q^(-4))e(t)
% 这里多项式A和C都只写出4阶，因为一般的经济时间序列阶数都不高。
clear
tic
% ====================第一步，模拟一个ARMA模型并绘制ACF，PACF图========================
%s首先设定ARMA模型的多项式系数。ARMA模型中只有多项式A(q)和C(q)，
%把A(q)的系数都设为0就得到MA模型，把C(q)的系数都设为0就得到AR模型。
a1 = -(0.5)^(1/3);
a2 = (0.5)^(2/3);
a3 = 0;
a4 = 0;
c1 = 0;
c2 = 0;
c3 = 0;
c4 = 0;
obv = 3000; %obv是模拟的观测数目。
A = [1 a1 a2 a3 a4];
B = []; %因为ARMA模型没有输入，因此多项式B是空的。
C = [1 c1 c2 c3 c4];
D = []; %把D也设为空的。
F = []; %ARMA模型里的F多项式也是空的。
m = idpoly(A,B,C,D,F,1,1) %这样就生成了ARMA模型，把它存储在m中。NoiseVariance被设定为1，1也是默认值。抽样间隔Ts设为1。
error = randn(obv,1); %生成一个obv*1的正态随机序列。准备用作模型的误差项。
e = iddata([],error,1); %用randn函数生成一个噪声序列。存储在e中。抽样间隔是1秒。
%u = []; %因为是ARMA模型，没有输出。所以把u设为空的。这句可省略。
y = sim(m,e);
get(y) %使用get函数来查看动态系统的所有性质。
r=y.OutputData; %把y.OutputData的全部值赋给变量r，r就是一个obv*1的向量。
figure(1)
plot(r) %绘出y随时间变化的曲线。或者写成plot(y)结果是一样的。
%==========================第二步，绘制ARMA序列r的ACF和PACF图=======================
%如果直接用函数autocorr和parcorr画图，总是会把滞后一阶的ACF和PACF也画出来，不好看，所以用下述方法画。
% figure(2)
% subplot(2,1,1)
% n=100;
% [ACF,Lags,Bounds]=autocorr(r,n,2);
% x=Lags(2:n);
% y=ACF(2:n); %注意这里的y和前面y的完全不同。
% h=stem(x,y,'fill','-');
% set(h(1),'Marker','.')
% hold on
% ylim([-1 1]);
% a=Bounds(1,1)*ones(1,n-1);
% line('XData',x,'YData',a,'Color','red','linestyle','--')
% line('XData',x,'YData',-a,'Color','red','linestyle','--')
% xlabel('lags')
% ylabel('ACF')
% title('ACF with 2stds Bounds(i.e., approximate 95% confidence interval)')
%
% subplot(2,1,2)
% [PACF,Lags,Bounds]=parcorr(r,n,2);
% x=Lags(2:n);
% y=PACF(2:n);
% h=stem(x,y,'fill','-');
% set(h(1),'Marker','.')
% hold on
% ylim([-1 1]);
% b=Bounds(1,1)*ones(1,n-1);
% line('XData',x,'YData',b,'Color','red','linestyle','--')
% line('XData',x,'YData',-b,'Color','red','linestyle','--')
% xlabel('lags')
% ylabel('PACF')
% title('PACF with 2stds Bounds(i.e., approximate 95% confidence interval)')
%%============================第三步，对r进行m阶时频归并========================
%%这一步对前面得到的ARMA序列r进行m阶时频归并。归并的方法是把每m个数据项加起来。
%%m的取值根据实际情况来决定。
m = 3; % m=3对应于把月度数据归并为季度数据。
R = reshape(r,m,obv/m); % 把向量 r 变形成m*(obv/m)的矩阵R.如果obv=3000,m=3，则R为3*1000的矩阵。
aggregatedr = sum(R); % sum(R)计算矩阵R每一列的和。得到的1*(obv/m)行向量aggregatedr就是时频归并后得到的序列。
dlmwrite('E:\Temporal Aggregation\Essays and Programs\Matlab code\output.txt',aggregatedr','delimiter','\t','precision',6,'newline','pc');
% 至此完成了对r的时频归并。
%%=====================第四步，画出aggregatedr的ACF和PACF图以判断其模型============
figure(3)
subplot(2,1,1)
n=100;
bound = 1;
[ACF,Lags,Bounds]=autocorr(aggregatedr,n,2);
x=Lags(2:n);
y=ACF(2:n); %注意这里的y和前面y的完全不同。
h=stem(x,y,'fill','-');
set(h(1),'Marker','.')
hold on
ylim([-bound bound]);
a=Bounds(1,1)*ones(1,n-1);
line('XData',x,'YData',a,'Color','red','linestyle','--')
line('XData',x,'YData',-a,'Color','red','linestyle','--')
xlabel('lags')
ylabel('ACF')
title('ACF with 2stds Bounds(i.e., approximate 95% confidence interval)')
subplot(2,1,2)
[PACF,Lags,Bounds]=parcorr(aggregatedr,n,2);
x=Lags(2:n);
y=PACF(2:n);
h=stem(x,y,'fill','-');
set(h(1),'Marker','.')
hold on
ylim([-bound bound]);
b=Bounds(1,1)*ones(1,n-1);
line('XData',x,'YData',b,'Color','red','linestyle','--')
line('XData',x,'YData',-b,'Color','red','linestyle','--')
xlabel('lags')
ylabel('PACF')
title('PACF with 2stds Bounds(i.e., approximate 95% confidence interval)')
t=toc
