Matlab應用：AR Model 預測台灣歷年人口數

自迴歸模型（Autoregressive model，簡稱AR模型），是統計上一種處理時間序列的方法，用同一變數例如 X 的之前各期，亦即 X1 至 Xt-1來預測本期 Xt的表現，並假設它們為一線性關係。因為這是從迴歸分析中的線性迴歸發展而來，只是不用 x預測 y，而是用 x 預測 x；所以叫做自迴歸。

其中： c 是常數項；被假設為平均數等於0且標準差等於的隨機誤差值；被假設為對於任何的 t 都不變。p 為 AR模型的階數。文字敘述為：X的當期值等於一個或數個落後期的線性組合，加常數項，加隨機誤差。將以上公式簡化為一階，方便解釋，也就是說我們利用 Xt-1 來預測 Xt ，其中 Xt 代表想要預測的當前時刻；Xt-1代表過去的時刻。而一階模型具有其特殊的意義存在，故我們給它一個專用的名字，稱為一階自我回歸模型(first-order autoregressive model)，簡稱為AR(1)，展開上式可得

 稱為一階自我回歸係數(first-order autoregressive coefficient)

clc
clear
close all

data_year = load('C:\TW_population_year.txt');
data_people = load('C:\TW_population_people.txt');
data_year_95 = data_year(1:61) % 提高難度，由95年預測105年
x = data_people(1:61)

x_avg = mean(x)

% 標準差
mystd = sqrt ( ( sum( (x - mean(x)).^2 ) )/ length(x) )
Mstd = std(x)

% 變異數
myvar =  ( sum( (x - mean(x)).^2 ) )/ length(x)
Mvar = var(x)

%% === AR Model 預測台灣歷年人口 ===
N = length(x)
for n=1:25
Y=x;
Y(1:n)=[];
m=N-n;
X=[];
for i=1:m
    for j=1:n
        X(i,j)=x(n+i-j);
    end
end
%待估計AR係數(利用最小二乘法求出ar model的係數，順便檢定)
fai=inv(X'*X)*(X'*Y)
%方差
Delta(n)=(Y-X*fai)'*(Y-X*fai)/(N-n);
%FPE準則
criterion(n,1)=(N+n)*Delta(1,n)/(N-n);
%AIC準則
criterion(n,2)=N*log(Delta(1,n))+2*n;
%BIC準則
criterion(n,3)=N*log(Delta(1,n))+n*log(log(N));
end

% 檢定結果
figure(1)
DELTA = plot(Delta')
hold on
FPE = plot(criterion(1:n,1))
legend('delta','FPE')

figure(2)
hold on
AIC = plot(criterion(1:n,2)) % 基本上看AIC與BIC，且以AIC為主，越低表示fitting越好
BIC = plot(criterion(1:n,3))
legend('AIC','BIC')

%% 運用AR Model預測(最小二乘法)
figure(3)
data = iddata(x,[]);
plot(data,'o')
hold on
arsys1 = [1;-fai(1:25)]'
arsys = idpoly([arsys1],[]);
predict_period = 10
y_pre = forecast(arsys,data,predict_period);
yh = y_pre.OutputData
yh = [data_people(1:61);yh]
plot(yh,'r'), legend('original','forecasted')
title('最小二乘法')

%% 預測(Levinson-Durbin遞推算法)
figure(4)
data = iddata(x,[]);
plot(data,'o')
hold on
arsys1 = aryule(x,25);
arsys = idpoly([arsys1],[]);
predict_period = 10
y_pre = forecast(arsys,data,predict_period);
yh = y_pre.OutputData;
yh = [data_people(1:61);yh]
plot(yh,'r'), legend('original','forecasted')
title('Levinson-Durbin遞推算法')

%% 預測(burg法)
figure(5)
data = iddata(x,[]);
plot(data,'o')
hold on
arsys1 = arburg(x,25);
arsys = idpoly([arsys1],[]);
predict_period = 10
y_pre = forecast(arsys,data,predict_period);
yh = y_pre.OutputData;
yh = [data_people(1:61);yh]
plot(yh,'r'), legend('original','forecasted')
title('Burg法')

由檢定結果可以看出，選用AR(25)應該會得到較好的fitting效果。比較了三種不同的AR建模法後發現，Burg法看似較為可信，遞推法基本上已經被我推到地下室去了 ..

沒想到Burg法預測十年後的人口數準確度可以到1%以內的誤差，不過這也可能是因為人口變動數據本身就不劇烈所以較好預測。 總之結果還算滿意! 重點是搞懂了該如何用 Matlab 實做AR Model了

補充說明：

有鑒於對於 AR Model 的數學式子沒有很明白的交待，在此我使用 Matlab help 中的 arburg 函數來輔助我們更了解 AR Model。 以下內容部分截錄至 help 文件中

Use a vector of polynomial coefficients to generate an AR(4) process by filtering 1024 samples of white noise. Reset the random number generator for reproducible results. Use Burg's method to estimate the coefficients.

淺見：建立一個 white noise 並且我們利用一個全極點的濾波器來進行濾波，再來利用 Burg's method 來逆推出這個全極點的濾波器參數為何。

剛開始看到這個範例時一直摸不著頭緒，究竟這樣的用處為何，後來我換個角度思考。假設我們手邊具有一個訊號，並且我們想利用 AR Model 來對這個訊號進行數學建模，但考量到訊號在量測時通常會參雜著 white noise ，以至於我們沒辦法觀測到原始訊號長什麼樣子，此時就可以利用Burg's method 對原始信號進行建模，消去 white noise。

以下為Matlab 裡面的範例

由此處我們可以發現，利用 filter 函數建立一個無零點，極點為 A 係數組成的轉移函式對 0.2*randn(1024,1) 進行濾波。再來利用 arburg 對 y 進行四階擬合，回傳的 arcoeffs 可以看出來與 A 基本上是相符的。

個人註記:
使用 AR model 進行 power spectrum 估計。

%%利用Yule--Walker方法進行功率譜估計
A = [1 -2.7607 3.8106 -2.6535 0.9238];

figure(7)
x = randn(1000,1);
y = filter(1,A,x);
[Pxx,F] = pyulear(y,4,1024,1); %% AR power spectrum 估計
%% [Pxx,w] = pyulear(y,4); plot(w*fs/2/pi,Pxx) % 若使用具有normalized
hold on
plot(F,10*log10(Pxx))
AR_par = aryule(data_people,4)

[Pxx,F] = pyulear(y,15,1024,1); %% AR power spectrum 估計
plot(F,10*log10(Pxx))
AR_par = aryule(data_people,15)

[H,F] = freqz(1,A,[],1);
plot(F,20*log10(abs(H)))

xlabel('Frequency (Hz)')
ylabel('PSD (dB/Hz)')
legend('pyulear PSD Estimate(4)','pyulear PSD Estimate(15)','True Power Spectral Density')

參考連結於Matlab應用：不專業的統計基礎札記下方