# Regression Analysis

### :small_blue_diamond: 回歸分析概述:
線性回歸是利用最小平方法，對自變數與應變數擬合出一條回歸線的一種統計方法。只有一個自變數是簡單回歸，當自變數有多項時則稱為多元回歸。如何在兩個變數當中找出一條直線，讓這條直接最接近變數之間變動的斜率，是19世紀統計學的主題。簡單回歸的方程式透過下列關係式來描述：<p class="text-center">
$y=\alpha x+\beta+\varepsilon,\varepsilon\sim N(0,\sigma^2)$

也就是當x增加α單位時，y會增加αx+β。換句話說只要知道這條回歸線的斜率與常數項，就能從x來預測y，而整個回歸分析的重點就在於找出α, β。其中<p class="text-center">
$\alpha=\frac{\sum_{j=1}^{n}{x_jy_j}- \sum_{j=1}^{n}{x_j}\sum_{j=1}^{n}{y_j}}{\sum_{j=1}^{n}{x_j^2}- (\sum_{j=1}^{n}{x_j})^2}$ ， $\beta=\frac{1}{n}(\sum_{j=1}^{n}{y_j - \alpha\sum_{j=1}^{n}{x_j}})$</p>

### :small_blue_diamond:簡單線性回歸分析:
Phillips曲線並不是直接以通貨膨脹率來分析，而是分析工資率漲跌與失業率之間的關係，因為在他的那篇經典論文中提到生活成本的調整會導致工資率的漲跌，進而影響勞動市場供需。我們仿照Phillips當年的作法，從中華民國統計資訊網彙整1990至2019年的工資變化率、失業率、物價指數成長率以及經濟成長率，來進行台灣版的菲利浦曲線。

讀取資料後，用對資料做線性回歸，將wage當作應變數、unemployment當作自變數建構工資率與失業率的回歸模型。從分析結果可看到，台版菲利浦曲線回歸模型為$y=−2.44x+11.75$。

接下來我們仿照Phillips當年的做法，畫一條台灣版的菲利浦曲線。從上述的模型可以發現這條回歸線的斜率為$-2.44$，斜率是負的代表工資變化率與失業率成反比：
#### 判定係數$R^2$
由於回歸模型是用最小平方法，擬合出一條回歸線，y能否被x完全解釋，取決於這條直線是不是能完全契合散布圖的形狀，也就是回歸模型的好壞。從散布圖中，我們可以知道回歸模型包含三種變異：總變異SST、可解釋變異SSR、其他變異SSE。聽起來有沒有一點熟悉？所以回歸模型的好壞，可以透過ANOVA來檢測。
既然已知總變異SST=可解釋變異SSR+其他變異SSE。那麼根據ANOVA就可計算模型的解釋力也就是判定係數$R^2$。<p class="text-center">
$R^2=\frac{SSR}{SST}=\frac{216.03}{261.03+147.76}=0.6385$</p>
從判定係數可以得知，回歸模型大約可以解釋64%的變異量，同時ANOVA F檢定顯著，回歸模型中的個別變數t檢定失業率也顯著，代表這個台灣版的菲利浦曲線確實呈現失業率與工資變化率的關係。(在簡單線性迴歸中，F test和t test的統計意義相同。判定係數: 迴歸模型的總變異中可被自變數解釋之百分比， 判定係數越大迴歸模型的配適度越好。一般而言，判定係數大於0.5就算不錯了。)
#### 圖形診斷
視覺診斷是最簡單判別殘差的方法，用plot()可以呼叫出4張殘差診斷圖。Residuals vs Fitted以及Scale-Location可以看出回歸模型與實際樣本的差異，Q-Q圖可以看出殘差是否符合常態分配，Residuals vs Leverage則可以透過「庫克距離」來衡量觀察值對模型的影響大小。
從4張殘差診斷圖中可以發現，台灣版菲利浦曲線的回歸模型，並沒有完全符合殘差常態分配的假定，殘差與樣本值的誤差也偏大，所以整體的模型解釋力只有0.6385。影響模型解釋力的主要因素，從4張診斷圖當中也可以發現那些樣本點是造成殘差的主要因素，也就是這些樣本屬於極端值。
#### 殘差檢定
和其他統計模型一樣，線性迴歸也有自己的前提條件:
(1)常態性(Normality) : 若母體資料呈現常態分配(Normal Distribution)，則誤差項也會呈現同樣的分配。可採用常態機率圖(normal probability plot) 或Shapiro-Wilk常態性檢定做檢查。
(2)獨立性(Independency) : 誤差項之間應該要相互獨立，否則在估計迴歸參數時會降低統計的檢定力。我們可以藉由Durbin-Watson test來檢查。
(3)變異數同質性(Constant Variance) : 變異數若不相等會導致自變數無法有效估計依變數。我們可以藉由殘差圖(Residual Plot)來檢查。

對以上範例做模型診斷，可以發現p值小於.05顯示不符合常態分配及殘差獨立性假設：

shapiro.test(model$residual) #殘差常態分配檢定
W = 0.80777, p-value = 0.0001155

durbinWatsonTest(model) #殘差常態自我相關(獨立性)檢定

| lag | Autocorrelation D-W | Statistic | p-value |
| --- | ------------------- | --------- |:------- |
| 1 | 0.2252316 | 1.227507 | 0.04 |

＊【小常識】
DW值愈接近2，代表其愈獨立。

ncvTest(model) #變異數同質性檢定
Non-constant Variance Score Test
Variance formula: ~ fitted.values
Chisquare = 0.3012812, Df = 1, p = 0.58308

### :small_blue_diamond: 多元線性回歸分析:
與前面自變數都只有一個的情況不同，多元線性回歸具有兩個以上的自變數，線性回歸方程式加入另一個變數之後成為：y=αx+βz+c。其實這已經不是一條直線，而是一個具有x、y、z軸的3D平面模型。我們現在將台灣版的菲利浦曲線，再加入一的自變數「經濟成長率」，來完成我們的多元回歸模型。

從上面分析顯示，多元回歸模型為y=−1.98x+0.40z+8.22，判定係數R^2=0.7。由回歸係數可得知失業率與工資變化率呈反比，也就是菲利浦曲線的情況，而經濟成長率與工資變化率成正比，相當符合直覺。引進scatterplot3d()後，我們可以描繪3個變數的3D立體關係。
#### 共線性
所有簡單線性回歸的模型診斷，包含殘差分析等，都適用多元回歸分析。這邊值得一提的是共線性(multicollinearity)問題，也就是兩個自變數是否高度相關。一般而言，當相關係數大於0.8時，代表變數關聯性強，這時只需要選擇其中一個納入回歸模型即可，否則會影響模型解釋力。

cor(phillips, use="complete.obs") #排除遺漏值後進行相關分析
| | year | wage | unemployment | cpi | growth |
| ------------ | ---------- | ---------- | ------------ | ---------- | ---------- |
| year | 1.0000000 | -0.6252511 | 0.7062332 | -0.5913482 | -0.5159854 |
| wage | -0.6252511 | 1.0000000 | -0.7990892 | 0.6952941 | 0.6246005 |
| unemployment | 0.7062332 | -0.7990892 | 1.0000000 | -0.7609087 | -0.5033919 |
| cpi | -0.5913482 | 0.6952941 | -0.7609087 | 1.0000000 | 0.4705541 |
| growth | -0.5159854 | 0.6246005 | -0.5033919 | 0.4705541 | 1.0000000 |

相關分析顯示unemployment與growth的相關係數是-0.503。
除此之外，變異數膨脹因子(Variance Inflation Factor, VIF)也是檢驗共線性的另一種方法。VIF應該越接近1越好，<5是可以接受的，但如果>10代表共線性的問題很大，不應該將變數納入模型，應該依據共線性大小逐一移除變數直到VIF<5為止。
模型VIF=1.339，共線性問題不大。
vif(MLR)
unemployment growth
1.339411 1.339411

| unemployment | growth |
| ------------ | -------- |
| 1.339411 | 1.339411 |

### :small_blue_diamond: 參考資料:
https://r-stat.neocities.org/regression
https://www.yongxi-stat.com/simple-regression-analysis/
https://www.yongxi-stat.com/multiple-regression-analysis/

# Time Series Analysis

### :small_blue_diamond: 什麼是時間序列分析:
時間序列分析（Time Series Analysis）是一種統計技術，用於分析隨時間變化的數據，旨在識別數據中的模式、趨勢和季節性，以便進行預測和決策。
### :small_blue_diamond: 時間序列的模式:
$Y_t=f(T_t,S_t,C_t,I_t)$

長期趨勢($T$)：指資料隨著時間遞移，呈現向上遞增或向下遞減的趨勢。例如物價指數的長期趨勢是呈現上揚。
循環變動($C$)：指資料呈現上下波動的情形。例如景氣循環有高峰、衰退、谷底、復甦的循環。
季節變動($S$)：指資料呈現依週、月、季規律變動的情形。例如旅遊業的旺季是寒暑假，淡季是冬天。
不規則變動($I$)：指資料隨機性的變動。不規則變動是除去長期趨勢、循環變動與季節變動後的隨機項。
依據上述四個因素，時間序列的數學表達式可分為相加模型(additive model)與相乘模型(multiplicative model)
1.相加模型:$Yt = Tt+St+Ct+IT$
2.相乘模型: $Yt = Tt\times St\times Ct\times IT$
相加模型假設各組成因素獨立互不影響，但在現實世界中任意因素很難相互獨立，所以經濟分析中多採用相乘模型。

### :small_blue_diamond: 分解方法(以乘法$Y_t = T_t \times S_t \times C_t \times I_t$為例):
1.運用移動平移法求出長期趨勢($T$)和週期變化($C$)，$MA=T\times C$
2.利用已知的$T\times C$ 推出 $\frac{Y}{T\times C} = S\times I$，接著因為發現 $I$ 可透過平均來消去，因此將$S\times I$做修正即可得 $S$ (季節性指數)
3.做圖選擇適合的曲線模型擬合序列的長期趨勢，即可得到 $T$ 長期趨勢
4.利用$\frac{T\times C}{T}$即可得 $C$ 週期變動因素
5.目前已的$S、T、C$ 即可推得 $I$ 不規則變動

### :small_blue_diamond: 時間序列預測:
預測是基於歷史時間資料在以後時間上進一步預測資料點的過程。
這可以使用統計模型來完成，例如:
1.自迴歸(AR)模型
2.移動平均(MA)模型
3.自迴歸移動平均(ARMA)模型
4.自迴歸綜合移動平均(ARIMA)模型

### :small_blue_diamond: 時間分析範例:
以遠東百貨為例，擷取2011年到2020年共計10年的營業收入進行分析。
原始資料:
time revenue(千元)
1 11-Jan 2493469
2 11-Feb 2215337
3 11-Mar 1400685
4 11-Apr 2393033
5 11-May 1878665
6 11-Jun 1737501

| Column 1 | time | revenue(千元) |
| -------- | ------ | ------------- |
| 1 | 11-Jan | 2493469 |
| 2 | 11-Feb | 2215337 |
| 3 | 11-Mar | 1400685 |
| 4 | 11-Apr | 2393033 |
| 5 | 11-May | 1878665 |
| 6 | 11-Jun | 1737501 |

此資料包含120筆記錄的時間序列資料。為進行時間序列分析，增加一個編號1到120的t時間變數，用來表示各時間的營業收入資料點：
time revenue t
1 11-Jan 2493469 1
2 11-Feb 2215337 2
3 11-Mar 1400685 3
4 11-Apr 2393033 4
5 11-May 1878665 5
6 11-Jun 1737501 6

|| time | revenue | t |
| -------- | ------ | ------- | --- |
| 1 | 11-Jan | 2493469 | 1 |
| 2 | 11-Feb | 2215337 | 2 |
| 3 | 11-Mar | 1400685 | 3 |
| 4 | 11-Apr | 2393033 | 4 |
| 5 | 11-May | 1878665 | 5 |
| 6 | 11-Jun | 1737501 | 6 |

繪製時間序列:

可以發現遠東百貨的營業收入有明顯季節起伏，長期趨勢在最近幾年些微下滑但還算平穩。

**長期趨勢($T$):**
要估計時間序列的長期趨勢，其實就是找出一條回歸線。這裡我們用簡單線性回歸來找出遠東百貨營業額的長期趨勢。
根據回歸模型，我們可以得出遠東百貨營業額長期趨勢線$y=1692t+3260076$。
從回歸模型當中可得知:隨著時間推移自2011年起遠東百貨每年營業額增加1,692千元。
**循環變動($C$):**
已知長期趨勢T後根據相乘模型$Y=T×C×S×I$，公式經過移項後循環變動可表示為：
$C×S×I=\frac{Y}{T}$
在沒有季節變動，例如年資料，或是不考慮季節變動的時候，我們可以刪除S，所以循環變動與不規模變動可表示為：
$C×I=\frac{Y}{T}$
循環變動($C$)與不規則變動($I$)是以實際營業額除以預測營業額來衡量。當數值大於1的時候，循環與不規則變動大於長期趨勢；反之則小於長期趨勢。在有季節變動的時候，循環與不規則變動可以與季節變動一起考量。因此接下來我們來分解季節變動。

季節變動
原始資料沒有季營收，只需將每3個月的營收相加就能計算出季營收。

[1]6109491 6009199 5101231 10420728 10266411 8772299 7754131 12335224
[9]11517708 10896653 10632521 13706045 11072973 10615980 10648762 13590122
[17] 11025307 10362100 10260305 13357659 10874476 10064518 9850916 12748750
[25] 9680792 9348923 9143721 12999668 9342711 8859550 8773507 12233554
[33] 9225858 8569486 8305841 11766113 8631254 7948294 8497408 12172036

計算出季營業額後，現在可以開始逐步拆解季節變動。根據相乘模型，季節變動與不規則變動可表示為：
$S×I=\frac{Y}{(T×C)}$
公式表明如果能計算出長期趨勢($T$)與循環變動($C$)，就能分離出季節與不規則變動。
那麼該如何計算長期趨勢與循環變動呢？一個簡單方法是計算季資料的移動平均，因為資料平滑後就相當於消除季節與不規則成分，只留下長期趨勢與循環變動。我們計算4期簡單移動平均數(利用四筆資料做平均接著往下滾一筆資料依此類推)可得:

[1] NA NA NA 6910162 7949392 8640167 9303392 9782016
[9] 10094841 10625929 11345527 11688232 11577048 11506880 11510940 11481959
[17] 11470043 11406573 11309459 11251343 11213635 11139240 11036892 10884665
[25] 10586244 10407345 10230547 10293276 10208756 10086413 9993859 9802331
[33] 9773117 9700601 9583685 9466825 9318174 9162876 9210767 9312248

一般而言為了資料對稱，移動平均期數都設定為奇數。這裡的期數4為偶數，所以接著再將兩兩移動平均相加，計算移動平均的移動平均，也就是中央移動平均數。最終可以獲得2011年第三季開始的季節中央移動平均數。這串移動平均數就是萃取出來的長期趨勢與循環變動。

[1] NA NA 7429777 8294780 8971780 9542704 9938428 10360385
[9] 10985728 11516879 11632640 11541964 11508910 11496450 11476001 11438308
[17] 11358016 11280401 11232489 11176437 11088066 10960779 10735455 10496795
[25] 10318946 10261911 10251016 10147584 10040136 9898095 9787724 9736859
[33] 9642143 9525255 9392499 9240525 9186821 9261508

比較原始季營收與移動平均調整後保留長期趨勢與循環變動的季營收如下：

接下來只要將季營業額除以中央移動平均數，就可換算出從2011年第三季開始到2020年第二季為止的季節變動與不規則變動百分比。

[1] NA NA 0.6865927 1.2562995 1.1443004 0.9192676 0.7802170
[8] 1.1906145 1.0484247 0.9461463 0.9140248 1.1874968 0.9621218 0.9234138
[15] 0.9279157 1.1881235 0.9707072 0.9185933 0.9134489 1.1951625 0.9807370
[22] 0.9182302 0.9176059 1.2145374 0.9381571 0.9110314 0.8919819 1.2810604
[29] 0.9305363 0.8950763 0.8963787 1.2564168 0.9568265 0.8996595 0.8843058
[36] 1.2733166 0.9395256 0.8582074 NA NA

現在已經可以歸納出遠東百貨每年Q1到Q4營業額的季節變動與不規則變動百分比：
Q1 Q2 Q3 Q4
2011 NA NA 0.6865927 1.256300
2012 1.1443004 0.9192676 0.7802170 1.190614
2013 1.0484247 0.9461463 0.9140248 1.187497
2014 0.9621218 0.9234138 0.9279157 1.188123
2015 0.9707072 0.9185933 0.9134489 1.195163
2016 0.9807370 0.9182302 0.9176059 1.214537
2017 0.9381571 0.9110314 0.8919819 1.281060
2018 0.9305363 0.8950763 0.8963787 1.256417
2019 0.9568265 0.8996595 0.8843058 1.273317
2020 0.9395256 0.8582074 NA NA

| | $Q_1$ | $Q_2$ | $Q_3$ | $Q_4$ |
|:---- |:--------- |:--------- |:--------- |:-------- |
| 2011 | NA | NA | 0.6865927 | 1.256300 |
| 2012 | 1.1443004 | 0.9192676 | 0.7802170 | 1.190614 |
| 2013 | 1.0484247 | 0.9461463 | 0.9140248 | 1.187497 |
| 2014 | 0.9621218 | 0.9234138 | 0.9279157 | 1.188123 |
| 2015 | 0.9707072 | 0.9185933 | 0.9134489 | 1.195163 |
| 2016 | 0.9807370 | 0.9182302 | 0.9176059 | 1.214537 |
| 2017 | 0.9381571 | 0.9110314 | 0.8919819 | 1.281060 |
| 2018 | 0.9305363 | 0.8950763 | 0.8963787 | 1.256417 |
| 2019 | 0.9568265 | 0.8996595 | 0.8843058 | 1.273317 |
| 2020 | 0.9395256 | 0.8582074 | NA | NA |

將每年Q1、Q2、Q3、Q4季節變動與不規則變動百分比平均，可獲得季節因子：

| Q1 | Q2 | Q3 | Q4 |
| --------- | --------- | --------- | --------- |
| 0.9857041 | 0.9099584 | 0.8680524 | 1.2270031 |

季節因子加總等於3.99，但四季總和其實應該等於4，所以我們做一下校正得出季節指數：
| Q1 | Q2 | Q3 | Q4 |
| --------- | --------- | --------- | --------- |
| 0.9879967 | 0.9120749 | 0.8700714| 1.2298570 |

到此季節變動已經拆解完成。下面這張表格就是最後的結果。我們可以看到第四季季節指數>1表示Q4是遠東百貨的旺季，而遠東百貨的淡季則落在Q3。

| | $Q_1$ | $Q_2$ | $Q_3$ | $Q_4$ |
|:-------- |:--------- |:--------- |:--------- |:--------- |
| 2011 | NA | NA | 0.6865927 | 1.256300 |
| 2012 | 1.1443004 | 0.9192676 | 0.7802170 | 1.190614 |
| 2013 | 1.0484247 | 0.9461463 | 0.9140248 | 1.187497 |
| 2014 | 0.9621218 | 0.9234138 | 0.9279157 | 1.188123 |
| 2015 | 0.9707072 | 0.9185933 | 0.9134489 | 1.195163 |
| 2016 | 0.9807370 | 0.9182302 | 0.9176059 | 1.214537 |
| 2017 | 0.9381571 | 0.9110314 | 0.8919819 | 1.281060 |
| 2018 | 0.9305363 | 0.8950763 | 0.8963787 | 1.256417 |
| 2019 | 0.9568265 | 0.8996595 | 0.8843058 | 1.273317 |
| 2020 | 0.9395256 | 0.8582074 | NA | NA |
| s_factor | 0.9857041 | 0.9099584 | 0.8680524 | 1.2270031 |
| s_index | 0.9879967 | 0.9120749 | 0.8700714 | 1.2298570 |

在經營管理策略上，有時候我們想瞭解營收成長單純是來自淡旺季的影響，還是每年都能維持穩定成長，這時候排除季節變動就顯得重要。有了季節指數就完成這項工作。只要將季營業額除以季節指數，就能得出排除季節與不規則變動影響的營業額。
[1] 6183716 6588493 5863003 8473122 10391139 9617959 8912063 10029803
[9] 11657638 11947103 12220286 11144422 11207500 11639373 12238952 11050164
[17] 11159255 11361019 11792486 10861148 11006591 11034750 11321963 10366042
[25] 9798405 10250170 10509162 10570065 9456217 9713621 10083663 9947135
[33] 9337944 9395595 9546160 9567058 8736116 8714519 9766334 9897115

經過季節校正後的營業額比較如下：

我們也可以建構時間序列，用R的ets()與auto.arima()模型自動預測遠東百貨未來的營業額。當然這個預測是不準確的，因為我們現在已經知道2021年受到covid-19衝擊各行各業營收下滑，這種黑天鵝事件是任何模型都無法預測的。

### :small_blue_diamond: 參考資料:
https://allaboutdataanalysis.medium.com/20%E5%80%8B%E6%99%82%E9%96%93%E5%BA%8F%E5%88%97%E5%9F%BA%E6%9C%AC%E6%A6%82%E5%BF%B5-%E5%AD%B8%E6%9C%83%E5%B9%AB%E4%BD%A0%E5%81%9A%E5%A5%BD%E9%A0%90%E6%B8%AC%E5%88%86%E6%9E%90-c297c631b374
https://r-stat.neocities.org/timeseries