Tugas STA1341 MPDW - Reviu Kasus Kaggle.com yang Menggunakan Analisis Regresi Peubah Lag
Pendahuluan
Tingkat pengangguran telah menjadi primary summary statistic untuk kesehatan pasar tenaga kerja selama beberapa waktu. Perkiraan tingkat pengangguran berperan saat pembuat kebijakan moneter mencoba merumuskan cara untuk mengkondisikan ekspektasi di lingkungan kebijakan baru.
Data yang digunakan
Data yang digunakan untuk analisis ini mewakili Statistik Pengangguran Area Lokal dari Januari 1990 - Desember 2016, dirinci menurut negara bagian dan bulan. Versi data dalam format CSV diperoleh dari Kaggle. Data mentah yang tidak diformat tersedia di Situs Web United States Bureau of Labor Statistics.
Catatan: Tingkat pengangguran ini adalah tarif U-3 bulanan umum dan tidak disesuaikan secara musiman atau dikategorikan berdasarkan usia, jenis kelamin, tingkat pendidikan, dll. (Penyesuaian musiman adalah teknik statistik yang mencoba mengukur dan menghilangkan pengaruh pola musiman yang dapat diprediksi untuk mengungkapkan bagaimana pekerjaan dan pengangguran berubah dari bulan ke bulan.)
Selanjutnya akan mengeksplorasi beberapa cara untuk menyesuaikan model yang dapat digunakan untuk meramalkan tingkat pengangguran di masa depan.
#working directory
setwd("E:/SEM 5/MPDW")
unemp.data <- read.csv("file.csv")
## Remove cases or rows with missing values. In this case we keep the rows which do not have NAs.
head(unemp.data[complete.cases(unemp.data), ])## Year Month State County Rate
## 1 2015 February Mississippi Newton County 6.1
## 2 2015 February Mississippi Panola County 9.4
## 3 2015 February Mississippi Monroe County 7.9
## 4 2015 February Mississippi Hinds County 6.1
## 5 2015 February Mississippi Kemper County 10.6
## 6 2015 February Mississippi Calhoun County 6.9tail(unemp.data[complete.cases(unemp.data), ])## Year Month State County Rate
## 885543 2009 November Maine Franklin County 10.2
## 885544 2009 November Maine Somerset County 10.5
## 885545 2009 November Maine Oxford County 10.5
## 885546 2009 November Maine Knox County 7.5
## 885547 2009 November Maine Piscataquis County 11.3
## 885548 2009 November Maine Aroostook County 9.0#Obtain general information about the dataset -
str(unemp.data)## 'data.frame': 885548 obs. of 5 variables:
## $ Year : int 2015 2015 2015 2015 2015 2015 2015 2015 2015 2015 ...
## $ Month : chr "February" "February" "February" "February" ...
## $ State : chr "Mississippi" "Mississippi" "Mississippi" "Mississippi" ...
## $ County: chr "Newton County" "Panola County" "Monroe County" "Hinds County" ...
## $ Rate : num 6.1 9.4 7.9 6.1 10.6 6.9 7.9 14.3 4.5 11.1 ...Terlihat bahwa dataset ini memiliki lebih dari 885 ribu pengamatan dari 5 peubah. Selanjutnya akan kita akan melihat pola yang diikuti oleh peubah rate dan melihat apakah dapat memodelkannya.
Plot rata-rata tingkat pengangguran tahunan nasional selama 27 tahun yang diteliti:
#Annual unemployment rate for the US from 1990 to 2016
Annual.rate <- sapply(split(unemp.data$Rate, unemp.data$Year), mean)
Year = c(1990:2016)
plot(Year, Annual.rate, main = "Annual mean US unemployment rate"); lines(Year, Annual.rate)
Dataset ini mencakup lebih dari 3 resesi. US mengalami resesi pada tahun 1990. Kemudian tak lama setelah George W. Bush menjabat, runtuhnya internet bubble dan serangan 11 September mendorong US mengalami resesi lain pada tahun 2001. Krisis keuangan 2007 hingga 2009 memunculkan resesi ketiga.
Time series Analysis
Analisis deret waktu merupakan metode untuk menganalisis data deret waktu. Peramalan deret waktu adalah penggunaan model untuk memprediksi nilai masa depan berdasarkan nilai yang diamati sebelumnya.
Mempersiapkan Dataset untuk analisis
Dataset ini memiliki hampir 900 ribu pengamatan dengan sekitar 3 ribu pengamatan merupakan data untuk setiap bulan dari tahun 1990 hingga 2016. Selanjutnya, akan membuat dataset lebih ringkas dengan membuat peubah tingkat pengangguran bulanan rata-rata untuk setiap 27 tahun.
#Break down into smaller datasets categorized by year
sub.1990 <- subset(unemp.data, unemp.data[,1] == 1990)
sub.1991 <- subset(unemp.data, unemp.data[,1] == 1991)
sub.1992 <- subset(unemp.data, unemp.data[,1] == 1992)
sub.1993 <- subset(unemp.data, unemp.data[,1] == 1993)
sub.1994 <- subset(unemp.data, unemp.data[,1] == 1994)
sub.1995 <- subset(unemp.data, unemp.data[,1] == 1995)
sub.1996 <- subset(unemp.data, unemp.data[,1] == 1996)
sub.1997 <- subset(unemp.data, unemp.data[,1] == 1997)
sub.1998 <- subset(unemp.data, unemp.data[,1] == 1998)
sub.1999 <- subset(unemp.data, unemp.data[,1] == 1999)
sub.2000 <- subset(unemp.data, unemp.data[,1] == 2000)
sub.2001 <- subset(unemp.data, unemp.data[,1] == 2001)
sub.2002 <- subset(unemp.data, unemp.data[,1] == 2002)
sub.2003 <- subset(unemp.data, unemp.data[,1] == 2003)
sub.2004 <- subset(unemp.data, unemp.data[,1] == 2004)
sub.2005 <- subset(unemp.data, unemp.data[,1] == 2005)
sub.2006 <- subset(unemp.data, unemp.data[,1] == 2006)
sub.2007 <- subset(unemp.data, unemp.data[,1] == 2007)
sub.2008 <- subset(unemp.data, unemp.data[,1] == 2008)
sub.2009 <- subset(unemp.data, unemp.data[,1] == 2009)
sub.2010 <- subset(unemp.data, unemp.data[,1] == 2010)
sub.2011 <- subset(unemp.data, unemp.data[,1] == 2011)
sub.2012 <- subset(unemp.data, unemp.data[,1] == 2012)
sub.2013 <- subset(unemp.data, unemp.data[,1] == 2013)
sub.2014 <- subset(unemp.data, unemp.data[,1] == 2014)
sub.2015 <- subset(unemp.data, unemp.data[,1] == 2015)
sub.2016 <- subset(unemp.data, unemp.data[,1] == 2016)#Prepare the data for time series analysis
#Mean Monthly employment rate sorted by year for the entire US
mean1990 <- sapply(split(sub.1990$Rate, factor(sub.1990$Month, month.name, ordered=TRUE)), mean)
mean1991 <- sapply(split(sub.1991$Rate, factor(sub.1991$Month, month.name, ordered=TRUE)), mean)
mean1992 <- sapply(split(sub.1992$Rate, factor(sub.1992$Month, month.name, ordered=TRUE)), mean)
mean1993 <- sapply(split(sub.1993$Rate, factor(sub.1993$Month, month.name, ordered=TRUE)), mean)
mean1994 <- sapply(split(sub.1994$Rate, factor(sub.1994$Month, month.name, ordered=TRUE)), mean)
mean1995 <- sapply(split(sub.1995$Rate, factor(sub.1995$Month, month.name, ordered=TRUE)), mean)
mean1996 <- sapply(split(sub.1996$Rate, factor(sub.1996$Month, month.name, ordered=TRUE)), mean)
mean1997 <- sapply(split(sub.1997$Rate, factor(sub.1997$Month, month.name, ordered=TRUE)), mean)
mean1998 <- sapply(split(sub.1998$Rate, factor(sub.1998$Month, month.name, ordered=TRUE)), mean)
mean1999 <- sapply(split(sub.1999$Rate, factor(sub.1999$Month, month.name, ordered=TRUE)), mean)
mean2000 <- sapply(split(sub.2000$Rate, factor(sub.2000$Month, month.name, ordered=TRUE)), mean)
mean2001 <- sapply(split(sub.2001$Rate, factor(sub.2001$Month, month.name, ordered=TRUE)), mean)
mean2002 <- sapply(split(sub.2002$Rate, factor(sub.2002$Month, month.name, ordered=TRUE)), mean)
mean2003 <- sapply(split(sub.2003$Rate, factor(sub.2003$Month, month.name, ordered=TRUE)), mean)
mean2004 <- sapply(split(sub.2004$Rate, factor(sub.2004$Month, month.name, ordered=TRUE)), mean)
mean2005 <- sapply(split(sub.2005$Rate, factor(sub.2005$Month, month.name, ordered=TRUE)), mean)
mean2006 <- sapply(split(sub.2006$Rate, factor(sub.2006$Month, month.name, ordered=TRUE)), mean)
mean2007 <- sapply(split(sub.2007$Rate, factor(sub.2007$Month, month.name, ordered=TRUE)), mean)
mean2008 <- sapply(split(sub.2008$Rate, factor(sub.2008$Month, month.name, ordered=TRUE)), mean)
mean2009 <- sapply(split(sub.2009$Rate, factor(sub.2009$Month, month.name, ordered=TRUE)), mean)
mean2010 <- sapply(split(sub.2010$Rate, factor(sub.2010$Month, month.name, ordered=TRUE)), mean)
mean2011 <- sapply(split(sub.2011$Rate, factor(sub.2011$Month, month.name, ordered=TRUE)), mean)
mean2012 <- sapply(split(sub.2012$Rate, factor(sub.2012$Month, month.name, ordered=TRUE)), mean)
mean2013 <- sapply(split(sub.2013$Rate, factor(sub.2013$Month, month.name, ordered=TRUE)), mean)
mean2014 <- sapply(split(sub.2014$Rate, factor(sub.2014$Month, month.name, ordered=TRUE)), mean)
mean2015 <- sapply(split(sub.2015$Rate, factor(sub.2015$Month, month.name, ordered=TRUE)), mean)
mean2016 <- sapply(split(sub.2016$Rate, factor(sub.2016$Month, month.name, ordered=TRUE)), mean)#Create a dataframe of mean monthly rates by year
mean.monthly.rate = data.frame('year'=c(rep(1990,12),rep(1991,12),rep(1992,12),rep(1993,12),rep(1994,12),rep(1995,12),rep(1996,12),rep(1997,12),rep(1998,12),rep(1999,12),rep(2000,12),rep(2001,12),rep(2002,12),rep(2003,12), rep(2004,12),rep(2005,12),rep(2006,12),rep(2007,12),rep(2008,12),rep(2009,12),rep(2010,12), rep(2011,12),rep(2012,12),rep(2013,12),rep(2014,12),rep(2015,12),rep(2016,12)),
'month' = rep(c('Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'),27),
'rate' = c(mean1990,mean1991,mean1992,mean1993,mean1994,mean1995,mean1996,mean1997,mean1998,mean1999, mean2000,mean2001,mean2002,mean2003,mean2004,mean2005,mean2006,mean2007,mean2008,mean2009,mean2010,mean2011,mean2012,mean2013,mean2014,mean2015,mean2016))#Look at the first few rows of this dataset
head(mean.monthly.rate)## year month rate
## 1 1990 Jan 7.055016
## 2 1990 Feb 7.123685
## 3 1990 Mar 6.516672
## 4 1990 Apr 6.006778
## 5 1990 May 5.614321
## 6 1990 Jun 5.704363str(mean.monthly.rate)## 'data.frame': 324 obs. of 3 variables:
## $ year : num 1990 1990 1990 1990 1990 1990 1990 1990 1990 1990 ...
## $ month: chr "Jan" "Feb" "Mar" "Apr" ...
## $ rate : num 7.06 7.12 6.52 6.01 5.61 ...Terlihat bahwa dengan menggabungkan data tiap bulan dari 27 tahun dengan menghitung rata-rata bulanannya dapat meringkas dataset yaitu mengurangi hampir 3 ribu pengamatan per bulan per tahun menjadi satu pengamatan per bulan per tahun dengan mengambil rata-rata bulanan. Terdapat 324 pengamatan dalam dataset baru ini dan selanjutnya akan digunakan untuk analisis lebih lanjut.
Pemodelan Regresi Deret Waktu
#Load required packages & functions
require(repr)
require(forecast)
#Convert to time series using the ts function
US.UnempRate.ts <- ts(mean.monthly.rate[,3], start = 1990, freq = 12)
#Plot the time series
options(repr.pmales.extlot.width=8, repr.plot.height=4)
plot(US.UnempRate.ts, main = 'Time Series Plot of US Unemployment Rate')
options(repr.pmales.extlot.width=8, repr.plot.height=4)
plot(log(US.UnempRate.ts), main = 'Time Series Plot of log US Unemployment Rate')
Data deret waktu ini memiliki tren yang signifikan, komponen musiman, dan juga memiliki pola siklus yang ditunjukkan oleh naik turunnya tingkat pengangguran yang sesuai dengan resesi dan ekspansi.Panjang siklus (dalam tahun) bervariasi dari waktu ke waktu. Plot dari deret log yang ditransformasi tampaknya memiliki besaran yang lebih merata untuk komponen musiman daripada yang tidak ditransformasi. Jadi akan menggunakan deret transformasi log untuk analisis lebih lanjut:
#Plot the distribution - histogram & QQ plot of the time series
dist.ts = function(df, col = 'residual', bins = 40)
{
par(mfrow = c(1,2))
temp = as.vector(df)
breaks = seq(min(temp), max(temp), length.out = (bins + 1))
hist(temp, breaks = breaks, main = paste('Distribution of', col), xlab = col)
qqnorm(temp, main = paste('Q-Q Normal plot of', col))
par(mfrow = c(1,1))
}
options(repr.pmales.extlot.width=8, repr.plot.height=4)
dist.ts(log(US.UnempRate.ts), col = 'US Unemp Rate')
Terlihat bahwa qq plot bukan garis lurus dan karenanya ini bukan i.i.d. distribusi normal, hal tersebut menegaskan temuan awal bahwa dataset ini memiliki pola dan tren tertentu.
Deteksi Autokorelasi
ACF adalah salah satu alat untuk menemukan pola dalam data, secara khusus dapat memberi tahu korelasi antara titik yang dipisahkan oleh berbagai jeda waktu. Untuk setiap deret, autokorelasi pada lag nol sama dengan satu. Selain itu, kita juga dapat mendefinisikan autokorelasi parsial orde kedua. Fungsi autokorelasi parsial (PACF) memberikan korelasi parsial deret waktu dengan nilai lag-nya sendiri, mengontrol nilai deret waktu di semua lag yang lebih pendek. Fungsi ini berperan penting dalam analisis data yang bertujuan untuk mengidentifikasi sejauh mana lag dalam model autoregressive (AR).
#Plot Auto correlation & Partial Auto Correlation funcions
plot.acf <- function(df, col = 'remainder', is.df =TRUE)
{
if(is.df) temp <- df[, col]
else temp <- df
par(mfrow = c(2,1))
acf(temp, main = paste('ACF of', col))
pacf(temp, main = paste('PACF of', col))
par(mfrow = c(1,1))
}
options(repr.pmales.extlot.width=8, repr.plot.height=6)
plot.acf(log(US.UnempRate.ts), col = 'US Unemp Rate', is.df = F)
Garis putus-putus biru pada plot ACF & PACF adalah interval kepercayaan 95% untuk nilai korelasi sedangkan garis vertikal menunjukkan korelasi. Setiap garis vertikal yang melintasi garis putus-putus berwarna biru menunjukkan nilai korelasi yang signifikan. Sebuah i.i.d. deret waktu yang terdistribusi normal tidak memiliki nilai ACF atau PACF yang signifikan di luar lag 0. Terlihat plot ACF di atas memiliki beberapa lag nonzero yang signifikan dan menunjukkan pola naik turun tertentu. Dari Plot PACF juga menunjukkan nilai yang signifikan untuk lebih dari satu lag, yang menunjukkan mungkin ada tren. Sehingga jelas ini bukan deret waktu stasioner.
Proses mengekstraksi komponen tren, siklus, dan musiman disebut sebagai dekomposisi. Selanjutnya, akan dilakukan metode dekomposisi.
Seasonal & Trend decomposition using Loess (STL Decomposition)
Dekomposisi Musiman dan Tren menggunakan model Loess atau STL, menguraikan deret waktu menjadi komponen musiman, tren, dan tidak teratur menggunakan loess dan memplot komponen secara terpisah, dimana komponen siklik termasuk dalam plot komponen tren. Selanjutnya, akan menerapkan model dekomposisi STL ke deret waktu transformasi log dan melihat apakah ini secara efektif menghilangkan tren dan musiman.
#Function for STL decomposition of a time series into components
ts.decomp <- function(df, col = 'Series Name', span = 0.13, Mult = TRUE)
{
if(Mult) temp = log(df) else temp = df
spans = span * length(temp)
fit <- stl(temp, s.window = "periodic", t.window = spans)
plot(fit, main = paste('Decompositon of',col,'with loess span = ', as.character(span)))
fit$time.series
}
US.Unemp.decomp <- ts.decomp(log(US.UnempRate.ts), col = 'US Unemp Rate', Mult = FALSE)
str(US.Unemp.decomp)## Time-Series [1:324, 1:3] from 1990 to 2017: 0.1571 0.1441 0.096 -0.0243 -0.0507 ...
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : chr [1:3] "seasonal" "trend" "remainder"#Plot the ACF & PACF of the remainder
options(repr.pmales.extlot.width=8, repr.plot.height=6)
plot.acf(US.Unemp.decomp[, 3], is.df = FALSE)
Dari plot ACF & PACF residual, masih terlihat beberapa nilai lag non-nol yang signifikan. Jadi proses STL belum sepenuhnya menghilangkan trend atau pola.
ARIMA Models
- Autoregressive Model
Nilai deret waktu autoregressive atau AR ditentukan oleh kombinasi linier dari nilai masa lalu. Model AR memperhitungkan korelasi serial dalam nilai deret waktu. Nilai deret autoregresif atau deret orde p atau AR(p) pada waktu t ditulis sebagai:
\(x_{t}=\alpha_{1}x_{t-1}+\alpha_{2}x_{t-2}+...+\alpha_{p}x_{t-p}+w_{t}\)
dengan :
ρ0=1
pk=αk
Jumlah nilai PACF bukan nol = p.
Model AR khusus untuk deret waktu stasioner, Jika varians tidak konstan atau komponen tren belum dihilangkan, model AR tidak akan menghasilkan hasil yang memuaskan.
- Moving Average Model
Untuk model rata-rata bergerak (MA), nilai deret waktu pada waktu t ditentukan oleh kombinasi linier dari istilah white noise masa lalu. Model MA memperhitungkan korelasi seri dalam istilah noise. \(x_{t}=w_{t}+\beta_{1}w_{t-1}+\beta_{2}w_{t-2}+...+\beta_{q}w_{t-q}\)
dengan:
ρ0=1
Number of nonzero ρk; nilai k≠0 = q.
Model MA juga khusus untuk deret waktu stasioner, jika varians tidak konstan model MA tidak akan menghasilkan hasil yang memuaskan. - Autoregressive Integrated Moving Average Model
Model ARIMA menambahkan istilah integrasi ke model ARMA. Komponen integrasi melakukan differencing untuk memodelkan komponen random walk. Komponen yang mengintegrasikan memodelkan salah satu bagian non-stasioner dari deret waktu. Model ARIMA didefinisikan oleh parameter p, d, q.
#ARIMA model estimation for STL remainder
ts.model = function(ts, col = 'remainder', order = c(0,0,1))
{
mod = arima(ts, order = order, include.mean = FALSE)
print(mod)
mod
}
arima.estimate1 <- ts.model(US.Unemp.decomp[, 3], order = c(1,0,1))#ARIMA(1,0,1) model##
## Call:
## arima(x = ts, order = order, include.mean = FALSE)
##
## Coefficients:
## ar1 ma1
## 0.8087 -0.0436
## s.e. 0.0382 0.0577
##
## sigma^2 estimated as 0.0007505: log likelihood = 705.32, aic = -1404.64arima.estimate2 <- ts.model(US.Unemp.decomp[, 3], order = c(0,1,1))#ARIMA(0,1,1) model##
## Call:
## arima(x = ts, order = order, include.mean = FALSE)
##
## Coefficients:
## ma1
## -0.1288
## s.e. 0.0516
##
## sigma^2 estimated as 0.0008241: log likelihood = 688.52, aic = -1373.05arima.estimate3 <- ts.model(US.Unemp.decomp[, 3], order = c(1,1,3))#ARIMA(1,1,3) model ##
## Call:
## arima(x = ts, order = order, include.mean = FALSE)
##
## Coefficients:
## ar1 ma1 ma2 ma3
## 0.7174 -0.9800 0.2555 -0.2754
## s.e. 0.0575 0.0701 0.0774 0.0696
##
## sigma^2 estimated as 0.0007203: log likelihood = 708.32, aic = -1406.63arima.estimate4 <- ts.model(US.Unemp.decomp[, 3], order = c(0,1,5))#ARIMA(0,1,5) model##
## Call:
## arima(x = ts, order = order, include.mean = FALSE)
##
## Coefficients:
## ma1 ma2 ma3 ma4 ma5
## -0.1625 -0.0053 -0.1550 -0.2964 -0.0006
## s.e. 0.0580 0.0614 0.0731 0.0618 0.0617
##
## sigma^2 estimated as 0.000757: log likelihood = 701.95, aic = -1391.89arima.estimate5 <- ts.model(US.Unemp.decomp[, 3], order = c(2,1,3))#ARIMA(2,1,3) model##
## Call:
## arima(x = ts, order = order, include.mean = FALSE)
##
## Coefficients:
## ar1 ar2 ma1 ma2 ma3
## 0.8605 -0.1226 -1.1136 0.3919 -0.2783
## s.e. 0.2490 0.2017 0.2403 0.2359 0.0687
##
## sigma^2 estimated as 0.0007191: log likelihood = 708.54, aic = -1405.09arima.estimate6 <- ts.model(US.Unemp.decomp[, 3], order = c(3,0,4))#ARIMA(3,0,4) model##
## Call:
## arima(x = ts, order = order, include.mean = FALSE)
##
## Coefficients:
## ar1 ar2 ar3 ma1 ma2 ma3 ma4
## -1.0076 0.2590 0.7093 1.8248 1.3592 0.5135 0.2456
## s.e. 0.0557 0.0958 0.0563 0.0696 0.1589 0.1603 0.0702
##
## sigma^2 estimated as 0.0006552: log likelihood = 725.9, aic = -1435.8arima.estimate7 <- ts.model(US.Unemp.decomp[, 3], order = c(4,0,5))#ARIMA(4,0,5) model##
## Call:
## arima(x = ts, order = order, include.mean = FALSE)
##
## Coefficients:
## ar1 ar2 ar3 ar4 ma1 ma2 ma3 ma4
## 1.0415 -0.2645 -0.1388 -0.0450 -0.2926 0.2524 0.1708 0.0666
## s.e. 0.1885 0.2460 0.2557 0.1586 0.1803 0.1862 0.1637 0.0779
## ma5
## 0.3459
## s.e. 0.0618
##
## sigma^2 estimated as 0.0006718: log likelihood = 722.95, aic = -1425.89arima.estimate8 <- ts.model(US.Unemp.decomp[, 3], order = c(4,1,5))#ARIMA(4,1,5) model##
## Call:
## arima(x = ts, order = order, include.mean = FALSE)
##
## Coefficients:
## ar1 ar2 ar3 ar4 ma1 ma2 ma3 ma4
## -0.8034 0.4454 0.6207 -0.1675 0.6331 -0.6143 -0.7050 -0.0755
## s.e. 0.3228 0.3061 0.1442 0.2595 0.3158 0.2735 0.2231 0.3362
## ma5
## -0.2382
## s.e. 0.0740
##
## sigma^2 estimated as 0.0006557: log likelihood = 722.13, aic = -1424.27cat(paste('Sigma^2 of the original series = ', as.character(var(log(US.UnempRate.ts)))))## Sigma^2 of the original series = 0.046737234933687Dari hasil di atas diperoleh bahwa Model ARIMA(3,0,4) dengan koefisien tersebut memiliki AIC (Akaike Information Criteria) paling kecil, Semua koefisien signifikan dan varians lebih kecil. Jadi model ini yang paling cocok untuk residual. Selanjutnya melihat plot ACF & PACF dan mengamati apakah komponen musiman & tren dari residual atau sisanya dihilangkan oleh model ini:
plot.acf(arima.estimate6$resid[-1], col = 'ARIMA(3,0,4) estimate', is.df = F)
Kedua plot memiliki nilai signifikan yang jauh lebih sedikit dibandingkan dengan sisa dekomposisi STL, menunjukkan ini mungkin merupakan perkiraan yang baik untuk STL remainder.
Forecasting
Paket R forecast memiliki alat untuk memperkirakan model ARIMA musiman yang optimal secara otomatis dan menggunakan model ini untuk membuat prakiraan di masa mendatang.
Fungsi auto.arima dari paket forecast melakukan regresi bertahap maju mulai dari urutan maksimum setiap koefisien untuk menemukan model ARIMA musiman yang optimal. Optimisasi dilakukan dengan menggunakan Akaike Information Criteria (AIC) dan Bayesian Information Criteria (BIC). Secara umum, model ARIMA dituliskan dengan notasi ARIMA (p d q), dimana p menyatakan orde dari proses autoregressive (AR), d menyatakan pembedaan (differencing),dan q menyatakan orde dari proses moving average (MA).
Selanjutnya akan menggunakan auto.arima untuk memperkirakan tingkat pengangguran di masa depan yaitu tingkat pengangguran untuk Januari hingga Desember 2017. Kami sudah memiliki data aktual untuk Januari hingga Juli 2017 dari situs web BLS, sehingga dapat melakukan perbandingkan rates yang diperkirakan dengan rates sebenarnya untuk melihat apakah model tersebut cocok atau tidak.
fit.US.Unemp <- auto.arima(log(US.UnempRate.ts), max.p = 5, max.q = 5, max.P = 2,
max.Q = 2, max.order = 5, max.d = 2, max.D = 2, start.p = 0,
start.q = 0, start.P = 0, start.Q = 0)
summary(fit.US.Unemp)## Series: log(US.UnempRate.ts)
## ARIMA(2,0,2)(1,1,1)[12]
##
## Coefficients:
## ar1 ar2 ma1 ma2 sar1 sma1
## 1.9448 -0.9486 -1.0487 0.1657 0.0484 -0.7496
## s.e. 0.0344 0.0338 0.0644 0.0575 0.0802 0.0540
##
## sigma^2 = 0.0008595: log likelihood = 655.65
## AIC=-1297.3 AICc=-1296.93 BIC=-1271.1
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.0002166452 0.02849136 0.02234327 -0.02780139 1.261069 0.2243726
## ACF1
## Training set -0.001567326Model yang diestimasi adalah ARIMA(2,0,2)(1,1,1). \(\sigma^{2}\) dari seri asli = 0,046737234933687 \(\sigma^{2}\) dari model yang diestimasi = 0,0008595 Standar error dari koefisien model umumnya urutan besarnya kurang dari nilai-nilai koefisien model. Semua koefisien model ini tampak signifikan.
unemp.forecast <- forecast(fit.US.Unemp, h=12)
summary(unemp.forecast)##
## Forecast method: ARIMA(2,0,2)(1,1,1)[12]
##
## Model Information:
## Series: log(US.UnempRate.ts)
## ARIMA(2,0,2)(1,1,1)[12]
##
## Coefficients:
## ar1 ar2 ma1 ma2 sar1 sma1
## 1.9448 -0.9486 -1.0487 0.1657 0.0484 -0.7496
## s.e. 0.0344 0.0338 0.0644 0.0575 0.0802 0.0540
##
## sigma^2 = 0.0008595: log likelihood = 655.65
## AIC=-1297.3 AICc=-1296.93 BIC=-1271.1
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.0002166452 0.02849136 0.02234327 -0.02780139 1.261069 0.2243726
## ACF1
## Training set -0.001567326
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2017 1.732220 1.694648 1.769791 1.674759 1.789681
## Feb 2017 1.721693 1.671245 1.772141 1.644539 1.798846
## Mar 2017 1.685626 1.623617 1.747634 1.590792 1.780459
## Apr 2017 1.571244 1.498421 1.644067 1.459871 1.682617
## May 2017 1.575594 1.492480 1.658709 1.448482 1.702707
## Jun 2017 1.647920 1.554939 1.740901 1.505717 1.790122
## Jul 2017 1.660378 1.557908 1.762847 1.503664 1.817092
## Aug 2017 1.621584 1.509985 1.733184 1.450907 1.792262
## Sep 2017 1.569445 1.449067 1.689823 1.385343 1.753548
## Oct 2017 1.556227 1.427423 1.685031 1.359238 1.753216
## Nov 2017 1.573184 1.436309 1.710060 1.363851 1.782518
## Dec 2017 1.635957 1.491369 1.780545 1.414829 1.857085plot(unemp.forecast)
Garis biru pada grafik di atas menunjukkan tingkat perkiraan dan warna abu-abu menunjukkan interval kepercayaan.Terlihat bahwa interval kepercayaan 80% dan 95% hampir sebanding dengan nilai perkiraan itu sendiri dan semakin lebar seiring berjalannya waktu, menggambarkan varians dapat meningkat seiring waktu, atau heteroskedastisitas. Itu berarti ini mungkin tidak cocok untuk data pengangguran.
Selanjutnya mari kita bandingkan nilai perkiraan dengan nilai aktual (tidak disesuaikan secara musiman) yang diperoleh dari situs web BLS secara berdampingan:
forecast.vs.actual = data.frame('year' = c(rep(2017,7)), 'month' = c('Jan','Feb','Mar','Apr','May','Jun','Jul'),
'forecasted'=c(round(exp(1.73222),2),round(exp(1.721693),2),round(exp(1.685626),2), round(exp(1.571244),2),round(exp(1.575594),2), round(exp(1.647920),2),round(exp(1.660378),2)),
'actual' = c(5.1,4.9,4.6,4.1,4.1,4.5,4.6))
head(forecast.vs.actual, 7)## year month forecasted actual
## 1 2017 Jan 5.65 5.1
## 2 2017 Feb 5.59 4.9
## 3 2017 Mar 5.40 4.6
## 4 2017 Apr 4.81 4.1
## 5 2017 May 4.83 4.1
## 6 2017 Jun 5.20 4.5
## 7 2017 Jul 5.26 4.6Ini menunjukkan bahwa seiring berjalannya waktu, nilai perkiraan semakin jauh dari aktual, menyiratkan model ini tidak cocok untuk data tersebut.
Kesimpulan
Tujuan analisis ini adalah untuk mempelajari tingkat pengangguran U-3 yang tidak disesuaikan secara musiman sebagai deret waktu univariat dan menemukan model linier yang akan menggambarkan variasi secara efektif. STL menguraikan deret waktu dan memasukkan model ARIMA(p,d,q) ke sisanya. Selain itu juga dilakukan penyesuaian model ARIMA(p,d,q)(P,D,Q) dan menggunakannya untuk meramalkan tingkat pengangguran untuk tahun 2017. Setelah membandingkan beberapa nilai aktual dengan nilai perkiraan, diperkirakan model ini kurang baik. Oleh karena itu, perlu analisis lebih lanjut atau dengan menyertakan beberapa faktor lain untuk memodelkan data tingkat pengangguran tersebut.
Sumber
Kamath Ranjani (2018). Time Series Analysis on US Unemployment Rate. Retrieved from https://www.kaggle.com/code/rkamath1/time-series-analysis-on-us-unemployment-rate/notebook