library(astsa)

# read data to R variable
birth.data<- read.csv("D:/Code/R/Time Serie/daily-total-female-births-in-cal.csv")
head(birth.data, 10)
# pull out number of births column
number_of_births<- birth.data$Daily.total.female.births.in.California..1959
number_of_births
##   [1] 35 32 30 31 44 29 45 43 38 27 38 33 55 47 45 37 50 43 41 52 34 53 39 32 37
##  [26] 43 39 35 44 38 24 23 31 44 38 50 38 51 31 31 51 36 45 51 34 52 47 45 46 39
##  [51] 48 37 35 52 42 45 39 37 30 35 28 45 34 36 50 44 39 32 39 45 43 39 31 27 30
##  [76] 42 46 41 36 45 46 43 38 34 35 56 36 32 50 41 39 41 47 34 36 33 35 38 38 34
## [101] 53 34 34 38 35 32 42 34 46 30 46 45 54 34 37 35 40 42 58 51 32 35 38 33 39
## [126] 47 38 52 30 34 40 35 42 41 42 38 24 34 43 36 55 41 45 41 37 43 39 33 43 40
## [151] 38 45 46 34 35 48 51 36 33 46 42 48 34 41 35 40 34 30 36 40 39 45 38 47 33
## [176] 30 42 43 41 41 59 43 45 38 37 45 42 57 46 51 41 47 26 35 44 41 42 36 45 45
## [201] 45 47 38 42 35 36 39 45 43 47 36 41 50 39 41 46 64 45 34 38 44 48 46 44 37
## [226] 39 44 45 33 44 38 46 46 40 39 44 48 50 41 42 51 41 44 38 68 40 42 51 44 45
## [251] 36 57 44 42 53 42 34 40 56 44 53 55 39 59 55 73 55 44 43 40 47 51 56 49 54
## [276] 56 47 44 43 42 45 50 48 43 40 59 41 42 51 49 45 43 42 38 47 38 36 42 35 28
## [301] 44 36 45 46 48 49 43 42 59 45 52 46 42 40 40 45 35 35 40 39 33 42 47 51 44
## [326] 40 57 49 45 49 51 46 44 52 45 32 46 41 34 33 36 49 43 43 34 39 35 52 47 52
## [351] 39 40 42 42 53 39 40 38 44 34 37 52 48 55 50
plot.ts(number_of_births ,main='Daily total female births in california, 1959', ylab = 'Number of births')

\[ \text{Box - Pierce test kiểm định giả thiết:} \]

\[ \begin{align*} H_0 &: \text{Không có sự tự tương quan trong chuỗi thời gian cho đến một độ trễ nhất định.} \\ &\rho(1)=\rho(2)=\ldots=\rho(h)=0 \\ \\ H_1 &: \text{Ít nhất có một sự tự tương quan trong chuỗi thời gian cho đến độ trễ được chỉ định khác 0.} \\ &\text{Ít nhất một } \rho(h) \text{ khác 0 với } h=1,2,\ldots,m \end{align*} \]

# Test for correlation
Box.test(number_of_births, lag = log(length(number_of_births)))
## 
##  Box-Pierce test
## 
## data:  number_of_births
## X-squared = 36.391, df = 5.8999, p-value = 2.088e-06
# Test for correlation in the differenced data
Box.test(diff(number_of_births), lag = log(length(diff(number_of_births))))
## 
##  Box-Pierce test
## 
## data:  diff(number_of_births)
## X-squared = 78.094, df = 5.8972, p-value = 7.661e-15
# Test for correlation in the differenced data
Box.test(diff(log(number_of_births)), lag = log(length(diff(number_of_births))))
## 
##  Box-Pierce test
## 
## data:  diff(log(number_of_births))
## X-squared = 72.32, df = 5.8972, p-value = 1.189e-13
diff_birth <-  diff(number_of_births)
diff_birth
##   [1]  -3  -2   1  13 -15  16  -2  -5 -11  11  -5  22  -8  -2  -8  13  -7  -2
##  [19]  11 -18  19 -14  -7   5   6  -4  -4   9  -6 -14  -1   8  13  -6  12 -12
##  [37]  13 -20   0  20 -15   9   6 -17  18  -5  -2   1  -7   9 -11  -2  17 -10
##  [55]   3  -6  -2  -7   5  -7  17 -11   2  14  -6  -5  -7   7   6  -2  -4  -8
##  [73]  -4   3  12   4  -5  -5   9   1  -3  -5  -4   1  21 -20  -4  18  -9  -2
##  [91]   2   6 -13   2  -3   2   3   0  -4  19 -19   0   4  -3  -3  10  -8  12
## [109] -16  16  -1   9 -20   3  -2   5   2  16  -7 -19   3   3  -5   6   8  -9
## [127]  14 -22   4   6  -5   7  -1   1  -4 -14  10   9  -7  19 -14   4  -4  -4
## [145]   6  -4  -6  10  -3  -2   7   1 -12   1  13   3 -15  -3  13  -4   6 -14
## [163]   7  -6   5  -6  -4   6   4  -1   6  -7   9 -14  -3  12   1  -2   0  18
## [181] -16   2  -7  -1   8  -3  15 -11   5 -10   6 -21   9   9  -3   1  -6   9
## [199]   0   0   2  -9   4  -7   1   3   6  -2   4 -11   5   9 -11   2   5  18
## [217] -19 -11   4   6   4  -2  -2  -7   2   5   1 -12  11  -6   8   0  -6  -1
## [235]   5   4   2  -9   1   9 -10   3  -6  30 -28   2   9  -7   1  -9  21 -13
## [253]  -2  11 -11  -8   6  16 -12   9   2 -16  20  -4  18 -18 -11  -1  -3   7
## [271]   4   5  -7   5   2  -9  -3  -1  -1   3   5  -2  -5  -3  19 -18   1   9
## [289]  -2  -4  -2  -1  -4   9  -9  -2   6  -7  -7  16  -8   9   1   2   1  -6
## [307]  -1  17 -14   7  -6  -4  -2   0   5 -10   0   5  -1  -6   9   5   4  -7
## [325]  -4  17  -8  -4   4   2  -5  -2   8  -7 -13  14  -5  -7  -1   3  13  -6
## [343]   0  -9   5  -4  17  -5   5 -13   1   2   0  11 -14   1  -2   6 -10   3
## [361]  15  -4   7  -5
# Plot the differenced data
plot.ts(diff(number_of_births), main='Differenced series', ylab = '')

\[ \text{Mô hình: } Y_t = p_t + Y_{t - 1} \\ Y_t \text{ là số lượng bé gái ra đời tại thời điểm t} \\ p_t \text{ là biến chênh lệch số lượng bé gái ra đời tại thời điểm t và } t - 1 \\ p_t = \delta(Y_t) = Y_t - Y_{t - 1} \]

# acf and pacf of the differenced data

acf(diff(number_of_births), main='ACF of differenced data', 50)

pacf(diff(number_of_births), main='PACF of differenced data', 50)

\[ \text{Xây dựng mô hình và ước lượng tham số} \\ \text{Mô hình MA(1): } \\ p_t = \theta_0 + \theta_1 \cdot u_{t - 1} + u_t \\ u_t \text{ là nhiễu trắng, thường là tuân theo phân bố chuẩn} \]

\[ \text{Mô hình MA(2): } \\ p_t = \theta_0 + \theta_1 \cdot u_{t - 1} + \theta_2 \cdot u_{t - 2} + u_t \]

\[ \text{Mô hình ARMA(7,1): } \\ p_t = \phi_1 \cdot p_{t - 1} + \phi_2 \cdot p_{t - 2} + \phi_3 \cdot p_{t - 3} + \phi_4 \cdot p_{t - 4} + \phi_5 \cdot p_{t - 5} + \phi_6 \cdot p_{t - 6} + \phi_7 \cdot p_{t - 7} + \theta_0 + \theta_1 \cdot u_{t - 1} + u_t \]

\[ \text{Mô hình ARMA(7,2): } \\ p_t = \phi_1 \cdot p_{t - 1} + \phi_2 \cdot p_{t - 2} + \phi_3 \cdot p_{t - 3} + \phi_4 \cdot p_{t - 4} + \phi_5 \cdot p_{t - 5} + \phi_6 \cdot p_{t - 6} + \phi_7 \cdot p_{t - 7} + \theta_0 + \theta_1 \cdot u_{t - 1} + \theta_2 \cdot u_{t - 2} + u_t \]

\[ \theta_0 \text{ là ước lượng giá trị trung bình của } p_t \]

model1<- arima(number_of_births, order=c(0,1,1))
SSE1<- sum(model1$residuals^2)
model1.test<- Box.test(model1$residuals, lag = log(length(model1$residuals)))

model2<- arima(number_of_births, order=c(0,1,2))
SSE2<- sum(model2$residuals^2)
model2.test<- Box.test(model2$residuals, lag = log(length(model2$residuals)))

model3<- arima(number_of_births, order=c(7,1,1))
SSE3<- sum(model3$residuals^2)
model3.test<- Box.test(model3$residuals, lag = log(length(model3$residuals)))

model4<- arima(number_of_births, order=c(7,1,2))
SSE4<- sum(model4$residuals^2)
model4.test<- Box.test(model4$residuals, lag = log(length(model4$residuals)))
df<- data.frame(row.names=c('AIC', 'SSE', 'p - value'), c(model1$aic, SSE1, model1.test$p.value), 
               c(model2$aic, SSE2, model2.test$p.value), c(model3$aic, SSE3, model3.test$p.value),
               c(model4$aic, SSE4, model4.test$p.value))
colnames(df)<- c('Arima(0,1,1)','Arima(0,1,2)', 'Arima(7,1,1)', 'Arima(7,1,2)')



format(df, scientific=FALSE)

\[ \text {Xét 4 mô hình trên, mô hình MA(2) có chỉ số AIC bé nhất.} \]

sarima(number_of_births, 0,1,2,0,0,0)
## initial  value 2.216721 
## iter   2 value 2.047518
## iter   3 value 1.974780
## iter   4 value 1.966955
## iter   5 value 1.958906
## iter   6 value 1.952299
## iter   7 value 1.951439
## iter   8 value 1.950801
## iter   9 value 1.950797
## iter  10 value 1.950650
## iter  11 value 1.950646
## iter  12 value 1.950638
## iter  13 value 1.950635
## iter  13 value 1.950635
## iter  13 value 1.950635
## final  value 1.950635 
## converged
## initial  value 1.950708 
## iter   2 value 1.950564
## iter   3 value 1.950290
## iter   4 value 1.950196
## iter   5 value 1.950185
## iter   6 value 1.950185
## iter   7 value 1.950185
## iter   7 value 1.950185
## iter   7 value 1.950185
## final  value 1.950185 
## converged

## $fit
## 
## Call:
## arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, Q), period = S), 
##     xreg = constant, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc, 
##         REPORT = 1, reltol = tol))
## 
## Coefficients:
##           ma1      ma2  constant
##       -0.8511  -0.1113     0.015
## s.e.   0.0496   0.0502     0.015
## 
## sigma^2 estimated as 49.08:  log likelihood = -1226.36,  aic = 2460.72
## 
## $degrees_of_freedom
## [1] 361
## 
## $ttable
##          Estimate     SE  t.value p.value
## ma1       -0.8511 0.0496 -17.1448  0.0000
## ma2       -0.1113 0.0502  -2.2164  0.0273
## constant   0.0150 0.0150   1.0007  0.3176
## 
## $AIC
## [1] 6.760225
## 
## $AICc
## [1] 6.760408
## 
## $BIC
## [1] 6.803051

\[ \text{Mô hình MA(2)} \\ p_t = \theta_0 + \theta_1 \cdot u_{t - 1} + \theta_2 \cdot u_{t - 2} + u_t \\ Y_t = Y_{t - 1} + \theta_0 + \theta_1 \cdot u_{t - 1} + \theta_2 \cdot u_{t - 2} + u_t \]

\[ \text{Ước lượng tham số: } \hat{\theta}_1 = -0.8511, \hat{\theta}_2 = -0.113, \hat{\theta}_0 = 0.015 \]

Mô hình ước lượng: \[ \hat{p_t} = 0.015 - 0.8511 \cdot u_{t - 1} - 0.113 \cdot u_{t - 2} + u_t \]

\[ \hat{Y_t} = Y_{t - 1} + 0.015 - 0.8511 \cdot u_{t - 1} - 0.113 \cdot u_{t - 2} + u_t \]

\[ \text{Ước lượng sai số chuẩn của các tham số: } SE(\hat{\theta}_1) = 0.0496, SE(\hat{\theta}_2) = 0.0502 \]

\[ \text{Ước lượng phương sai của } u_t: \sigma^2 = 49.08 \]

\[ \text{Kiểm định các hệ số } \theta_1, \theta_2, \theta_0 \text { có thực sự khác khác 0} \\ H_0: \theta_1 = 0 \\ H_1: \theta_1 \text{ khác 0} \\ \]

\[ \text{Do } p_{\text{value}} < 0.05 \text{ nên bác bỏ } H_0. \]

\[ \text{ Có thể nói } \theta_1 \text{ thực sự khác 0.} \]

\[ H_0: \theta_2 = 0 \\ H_1: \theta_2 \text{ khác 0} \\ \]

\[ \text{Do } p_{\text{value}} < 0.05 \text{ nên bác bỏ } H_0. \]

\[ \text{ Có thể nói } \theta_2 \text{ thực sự khác 0.} \]

\[ H_0: \theta_0 = 0 \\ H_1: \theta_0 \text{ khác 0} \\ \]

\[ \text{Do } p_{\text{value}} > 0.05 \text{ nên không có cơ sở bác bỏ } H_0. \]

\[ \text{Ljung - Box statistic:} \\ \text{Kiểm định xem phần dư có được phân phối độc lập không?} \\ H_0: \text{Phần dư được phân phối độc lập} \\ H_1: \text{Phần dư không được phân phối độc lập, tức là chúng có tương quan.} \\ p_{\text{value}} > 0.05. \text{ Chấp nhận } H_0, \text{phần dư được phân phối độc lập.} \]

# Load necessary libraries
library(forecast)
library(ggplot2)

# Update the model with the estimated sigma^2
model2$sigma2 <-  49.08

# Forecast the next 30 days
forecast_values <-  forecast(model2, h=60)

# Plot the original data and the forecast
autoplot(forecast_values) + 
  ggtitle("Daily total female births in California, 1959") +
  xlab("Time") + 
  ylab("Number of births")