Gọi packages

library(readxl)
library(fBasics)
## Loading required package: timeDate
## Loading required package: timeSeries
library(tseries)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(ggplot2)
library(rmarkdown)

Đọc dữ liệu

Để đánh giá hiệu quả kinh doanh của công ty … giai đoạn …, cần đánh giá được thực trạng tình hình kinh doanh của công ty qua 3 loại số liệu chính: Doanh thu – Chi phí và Vốn kinh doanh.

Thu thập số liệu qua báo cáo Kinh doanh của công ty, ta thu được bảng dữ liệu được tổng hợp trong giai đoạn … như sau:

Bảng 1: Doanh thu - chi phí và vốn kinh doanh của Công ty … giai đoạn …

dat <- read_excel("Số-liệu-công-ty-CP-6-NA.xlsx")
paged_table(dat)
rev <- dat$doanh_thu
time <- dat$time
lrev <- log(rev) # lấy logarithm

Dữ liệu thu thập được tại Bảng 1 là dạng dữ liệu được thống kê theo quý về doanh thu của Công ty … trong khoảng thời gian từ quý …

Có thể thấy được sự tăng trưởng vượt bậc của công ty trong kỳ nghiên cứu, từ quý …, doanh thu đã tăng gấp 9 lần. Biến động về tổng doanh thu của Công ty TNHH Truyền thông Thịnh Vượng theo từng quý thông qua biểu đồ doanh thu:

library("zoo")
## 
## Attaching package: 'zoo'
## The following object is masked from 'package:timeSeries':
## 
##     time<-
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
time <- as.yearqtr(time, format = "Q%q.%Y")

ggplot(NULL) + 
  aes(x = time, y = rev/1000000000) + 
  geom_line(col = "124B52") + 
  ylab("Doanh thu (tỷ đồng)") +
  xlab("Quý") +
  theme_classic()

Tính từ quý 1 đến hết giai đoạn, doanh thu tăng hơn 6 lần

Giải thích vì sao doanh thu tăng trưởng lớn trong giai đoạn này (lưu ý vấn đề covid?)

Ta có đồ thị ở bên dưới thể hiện diễn biến doanh thu đã được lấy logarit tự nhiên (việc lấy loga nepe để giảm sự biến thiên của biến doanh thu, nhưng vẫn giữ được xu hướng của chuỗi) theo thời gian của Công ty … giai đoạn …

ggplot(NULL) + 
  aes(x = time, y = lrev) + 
  geom_line(col = "124B52") + 
  ylab("log(Doanh thu)") +
  xlab("Quý") +
  theme_classic()

Ta thấy chuỗi giá trị có xu thế tăng dần đều (có vẻ không dừng), tuy nhiên, để xác định chuỗi có đảm bảo tính dừng hay không, cần kiểm định Unitroot để xác định

Kiểm định ADF cho chuỗi:

library(tseries)
adf.test(lrev)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  lrev
## Dickey-Fuller = -0.71786, Lag order = 2, p-value = 0.9568
## alternative hypothesis: stationary

Kiểm định cho p-value > 5% => chưa đủ cơ sở bác bỏ H0 (có nghiệm đơn vị) tức chuỗi không dừng

Để xử lý vấn đề chuỗi không dừng, lấy sai phân của chuỗi (tính ∆y) để chuyển chuỗi không dừng (y) sang chuỗi dừng (∆y)

Đồ thị ở bên hình dưới ta thấy nó đã biến động ổn định xung quanh đường trung bình (có thể chuỗi này đã dừng)

dlrev <- diff(lrev)
ggplot(NULL) + 
  aes(x = time[-1], y = dlrev) + 
  geom_line(col = "124B52") + 
  ylab("∆log(Doanh thu)") +
  xlab("Quý") +
  theme_classic()

Kiểm tra tính dừng của chuỗi đã lấy sai phân (∆log(Doanh thu)):

adf.test(dlrev)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  dlrev
## Dickey-Fuller = -4.0097, Lag order = 2, p-value = 0.02292
## alternative hypothesis: stationary

Ở mức 5%: ta thấy p-value = 0.022 -> bác bỏ H0 -> chuỗi đã dừng đã có thể áp dụng phương pháp AR, MA

ARIMA

Kiểm tra đồ thị ACF và PACF

acf(dlrev)

pacf(dlrev)

  1. Nhìn quan đồ thị tự tương quan ACF: ta thấy chuỗi có tự tương quan cao ở các bậc trong quá khứ (bậc 0 cao - em nghĩ là tương quan cao với y ở hiện tại)

  2. Nhìn quan đồ thị PACF: ta thấy chuỗi có tự tương quan riêng cao ở bậc 2

Từ 1 và 2. Đề xuất chọn mô hình AR(2) tức mô hình AR với bậc 2 để phân tích

Chạy mô hình ARIMA(2, 1, 0) cho log(doanh thu) hay AR(2) cho ∆log(doanh thu)

m <- arima(dlrev, order = c(2, 0, 0), method = "ML")
m
## 
## Call:
## arima(x = dlrev, order = c(2, 0, 0), method = "ML")
## 
## Coefficients:
##           ar1      ar2  intercept
##       -0.3616  -0.5211     0.1593
## s.e.   0.2337   0.2182     0.0361
## 
## sigma^2 estimated as 0.05748:  log likelihood = -0.22,  aic = 8.43

Thấy được hệ số của ar2 có tác động âm, hệ số aic đạt 8.43

Tính p-value cho các hệ số trong mô hình

format((1-pnorm(abs(m$coef)/sqrt(diag(m$var.coef))))*2, scientific = F)
##             ar1             ar2       intercept 
## "0.12180762815" "0.01692012874" "0.00001000416"

Hệ số ar2 có p-value = 0,0169 như vậy đạt ý nghĩa mức 5%

Kiểm định mô hình ARIMA

Kiểm tra phần dư của mô hình đã nhiễu trắng hay chưa

Một chuỗi nhiễu trắng là chuỗi có trung bình = 0, phương sai không đổi, và không có hiện tương tự tương quan

=> Tức chuỗi nhiễu trắng là chuỗi dừng xung quanh giá trị trung bình = 0, và không có hiện tương tự tương quan

Kiểm tra tự tương quan phần dư

plot(m$residuals)

Nhìn đồ thị phần dư bên trên ta thấy có vẻ đã dừng xung quanh đường trung bình là 0

Đồ thị bên dưới thể hiện phần dư có phân phối chuẩn hay không (đây cũng là 1 điều kiện trong mô hình)

Cụ thể, nếu các điểm nằm thành 1 đường thẳng từ dưới lên => phần dư theo phân phối chuẩn

qqplot(m$residuals, dlrev)

Kiểm định tính dừng của phân dư mô hình bằng kiểm định ADF (Augmented Dickey-Fuller Test)

adf.test(m$residuals)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  m$residuals
## Dickey-Fuller = -4.0671, Lag order = 2, p-value = 0.02092
## alternative hypothesis: stationary

p-value < 5% => Đủ cơ sở để bác bỏ H0, tức phần dư của mô hình đã dừng

Tiếp theo ta kiểm định phần dư của mô hình có phân phối chuẩn hay không bằng kiểm định Jarque - Bera

normalTest(m$residuals, method = "jb")
## 
## Title:
##  Jarque - Bera Normalality Test
## 
## Test Results:
##   STATISTIC:
##     X-squared: 0.5839
##   P VALUE:
##     Asymptotic p Value: 0.7468 
## 
## Description:
##  Tue Mar 29 22:04:48 2022 by user:

Kiểm định tự tương quan cho phân dư bằng kiểm định Box - Pierce tới bậc trễ thứ 10 (em lấy đại con số 10 này - chỉ cần lấy đủ xa với hiện tại là được)

Box.test(m$residuals, lag = 10)
## 
##  Box-Pierce test
## 
## data:  m$residuals
## X-squared = 2.5757, df = 10, p-value = 0.9897

Từ 2 đồ thị về phần dư, Kiểm định ADF và Kiểm định Box - Pierce (p-value > 5% => chưa đủ cơ sở bác bỏ H0 => không có tự tương quan)

=> Phần dư đã dừng và không có tự tương quan => nhiễu trắng

=> Lưu ý: vì trong trường hợp này phần dư không có hiện tượng tự tương quan nữa -> không cần tăng bậc q cho mô hình bình quân trượt MA

Ngoài ta kiểm định Jarque - Bera còn cho thấy rằng phần dư đã tuân theo quy luận phân phối chuẩn (thỏa điều kiện)

Dự báo cho 5 quý tiếp theo

pred <- predict(m, 5)
pred
## $pred
## Time Series:
## Start = 15 
## End = 19 
## Frequency = 1 
## [1] 0.35869383 0.18958699 0.04441771 0.18503076 0.20982692
## 
## $se
## Time Series:
## Start = 15 
## End = 19 
## Frequency = 1 
## [1] 0.2397549 0.2549494 0.2715808 0.2828422 0.2835618
`Doanh thu dự báo` <- c(
  exp(pred$pred[1] + lrev[15]),
  exp(pred$pred[2] + lrev[15] + pred$pred[1]),
  exp(pred$pred[3] + lrev[15] + pred$pred[2]),
  exp(pred$pred[4] + lrev[15] + pred$pred[3]),
  exp(pred$pred[5] + lrev[15] + pred$pred[4])
)
time <- c(
  "Q4.2021",
  "Q1.2022",
  "Q2.2022",
  "Q3.2022",
  "Q4.2022"
)
pred_rev <- data.frame(time, `Doanh thu dự báo`)
paged_table(pred_rev)

Phân tích xu hướng biến động của nhóm chỉ tiêu hiệu quả kinh doanh.

chitieu <- read_excel("Số-liệu-công-ty-CP-6-NA.xlsx", sheet = 2)
paged_table(chitieu)

Đồ thị gốc

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:timeSeries':
## 
##     filter, lag
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)
library(zoo)
chitieu %>% gather(
    doanhthu_chiphi,
    doanhthu_vonkd,
    doanhloi_chiphi, 
    doanhloi_doanhthu,
    doanhloi_vonkd,
    key = ten, 
    value = value
  ) %>%
  mutate(ten = case_when(
    ten == "doanhthu_chiphi" ~ "Doanh thu trên 1 đồng chi phí sản xuất",
    ten == "doanhthu_vonkd" ~ "Doanh thu trên 1 đồng vốn kinh doanh",
    ten == "doanhloi_chiphi" ~ "Doanh lợi theo chi phí sản xuất",
    ten == "doanhloi_doanhthu" ~ "Doanh lợi theo doanh thu",
    ten == "doanhloi_vonkd" ~ "Doanh lợi theo vốn kinh doanh"
  )) %>%
  arrange(time, desc(ten)) %>%
  mutate(time = as.yearqtr(time, format = "Q%q.%Y")) %>%
    ggplot(aes(time, value, col = ten)) +
    geom_line() +
    guides(col=guide_legend("Chỉ tiêu tài chính")) +
    theme_classic()+
    ylab("") +
    xlab("")

Đồ thị san bằng mũ với alpha = 0.25

library(fpp2)
## ── Attaching packages ────────────────────────────────────────────── fpp2 2.4 ──
## ✓ forecast  8.16     ✓ expsmooth 2.3 
## ✓ fma       2.4
## 
ses1 <- ses(chitieu$doanhthu_chiphi, alpha = 0.25)
ses2 <- ses(chitieu$doanhthu_vonkd, alpha = 0.25)
ses3 <- ses(chitieu$doanhloi_chiphi, alpha = 0.25)
ses4 <- ses(chitieu$doanhloi_vonkd, alpha = 0.25)
ses5 <- ses(chitieu$doanhloi_doanhthu, alpha = 0.25)

chitieu %>% mutate(
    doanhthu_chiphi = ses1$fitted,
    doanhthu_vonkd = ses2$fitted,
    doanhloi_chiphi = ses3$fitted,
    doanhloi_vonkd = ses4$fitted,
    doanhloi_doanhthu = ses5$fitted
  ) %>% gather(
    doanhthu_chiphi, 
    doanhthu_vonkd, 
    doanhloi_chiphi, 
    doanhloi_doanhthu, 
    doanhloi_vonkd, 
    key = ten, 
    value = value
  ) %>%
  mutate(ten = case_when(
    ten == "doanhthu_chiphi" ~ "Doanh thu trên 1 đồng chi phí sản xuất",
    ten == "doanhthu_vonkd" ~ "Doanh thu trên 1 đồng vốn kinh doanh",
    ten == "doanhloi_chiphi" ~ "Doanh lợi theo chi phí sản xuất",
    ten == "doanhloi_doanhthu" ~ "Doanh lợi theo doanh thu",
    ten == "doanhloi_vonkd" ~ "Doanh lợi theo vốn kinh doanh"
  )) %>%
  arrange(time, desc(ten)) %>%
  mutate(time = as.yearqtr(time, format = "Q%q.%Y")) %>%
  ggplot(aes(time, value, col = ten)) +
  geom_line() + 
  guides(col=guide_legend("Chỉ tiêu tài chính")) +
  theme_classic()+
  ylab("") + 
  xlab("")

Đồ thị san bằng mũ với alpha tốt nhất do máy tính tự tính toán

ses1 <- ses(chitieu$doanhthu_chiphi)
ses2 <- ses(chitieu$doanhthu_vonkd)
ses3 <- ses(chitieu$doanhloi_chiphi)
ses4 <- ses(chitieu$doanhloi_vonkd)
ses5 <- ses(chitieu$doanhloi_doanhthu)

chitieu %>% mutate(
    doanhthu_chiphi = ses1$fitted,
    doanhthu_vonkd = ses2$fitted,
    doanhloi_chiphi = ses3$fitted,
    doanhloi_vonkd = ses4$fitted,
    doanhloi_doanhthu = ses5$fitted
  ) %>% gather(
    doanhthu_chiphi, 
    doanhthu_vonkd, 
    doanhloi_chiphi, 
    doanhloi_doanhthu, 
    doanhloi_vonkd, 
    key = ten, 
    value = value
  ) %>%
  mutate(ten = case_when(
    ten == "doanhthu_chiphi" ~ "Doanh thu trên 1 đồng chi phí sản xuất",
    ten == "doanhthu_vonkd" ~ "Doanh thu trên 1 đồng vốn kinh doanh",
    ten == "doanhloi_chiphi" ~ "Doanh lợi theo chi phí sản xuất",
    ten == "doanhloi_doanhthu" ~ "Doanh lợi theo doanh thu",
    ten == "doanhloi_vonkd" ~ "Doanh lợi theo vốn kinh doanh"
  )) %>%
  arrange(time, desc(ten)) %>%
  mutate(time = as.yearqtr(time, format = "Q%q.%Y")) %>%
  ggplot(aes(time, value, col = ten)) +
  geom_line() + 
  guides(col=guide_legend("Chỉ tiêu tài chính")) +
  theme_classic()+
  ylab("") + 
  xlab("")

Dự đoán đến hết năm 2022 bằng phương pháp Holt - Winter có xu thế dạng cộng

Doanh thu trên 1 đồng chi phí sản xuất

hwa1 <- hw(
  ts(chitieu$doanhthu_chiphi, frequency = 4, start = c(2018, 1), end = c(2021, 3)), 
  h = 5, 
  seasonal = "additive",
  level = 95
)
paged_table(data.frame(hwa1))

Doanh thu trên 1 đồng vốn kinh doanh

hwa2 <- hw(
  ts(chitieu$doanhthu_vonkd, frequency = 4, start = c(2018, 1), end = c(2021, 3)), 
  h = 5, 
  seasonal = "additive",
  level = 95
)
paged_table(data.frame(hwa2))

Doanh lợi theo chi phí sản xuất

hwa3 <- hw(
  ts(chitieu$doanhloi_chiphi, frequency = 4, start = c(2018, 1), end = c(2021, 3)), 
  h = 5, 
  seasonal = "additive",
  level = 95
)

paged_table(data.frame(hwa3))

Doanh lợi theo doanh thu

hwa4 <- hw(
  ts(chitieu$doanhloi_doanhthu, frequency = 4, start = c(2018, 1), end = c(2021, 3)), 
  h = 5, 
  seasonal = "additive",
  level = 95
)

paged_table(data.frame(hwa4))

Doanh lợi theo vốn kinh doanh

hwa5 <- hw(
  ts(chitieu$doanhloi_vonkd, frequency = 4, start = c(2018, 1), end = c(2021, 3)), 
  h = 5, 
  seasonal = "additive",
  level = 95
)
paged_table(data.frame(hwa5))

Dự đoán đến hết năm 2022 bằng phương pháp Holt - Winter có xu thế dạng nhân

Doanh thu trên 1 đồng chi phí sản xuất

hwm1 <- hw(
  ts(chitieu$doanhthu_chiphi, frequency = 4, start = c(2018, 1), end = c(2021, 3)), 
  h = 5, 
  seasonal = "multiplicative",
  level = 95
)

paged_table(data.frame(hwm1))

Doanh thu trên 1 đồng vốn kinh doanh

hwm2 <- hw(
  ts(chitieu$doanhthu_vonkd, frequency = 4, start = c(2018, 1), end = c(2021, 3)), 
  h = 5, 
  seasonal = "multiplicative",
  level = 95
)

paged_table(data.frame(hwm2))

Doanh lợi theo chi phí sản xuất

hwm3 <- hw(
  ts(chitieu$doanhloi_chiphi, frequency = 4, start = c(2018, 1), end = c(2021, 3)), 
  h = 5, 
  seasonal = "multiplicative",
  level = 95
)

paged_table(data.frame(hwm3))

Doanh lợi theo doanh thu

hwm4 <- hw(
  ts(chitieu$doanhloi_doanhthu, frequency = 4, start = c(2018, 1), end = c(2021, 3)), 
  h = 5, 
  seasonal = "multiplicative",
  level = 95
)

paged_table(data.frame(hwm4))

Doanh lợi theo vốn kinh doanh

hwm5 <- hw(
  ts(chitieu$doanhloi_vonkd, frequency = 4, start = c(2018, 1), end = c(2021, 3)), 
  h = 5, 
  seasonal = "multiplicative",
  level = 95
)

paged_table(data.frame(hwm5))
RMSE <- data.frame(`RMSE xu thế dạng cộng` = c(
  accuracy(hwa1)[, 2],
  accuracy(hwa2)[, 2],
  accuracy(hwa3)[, 2],
  accuracy(hwa4)[, 2],
  accuracy(hwa5)[, 2]
), 
`RMSE xu thế dạng nhân` = c(
  accuracy(hwm1)[, 2],
  accuracy(hwm2)[, 2],
  accuracy(hwm3)[, 2],
  accuracy(hwm4)[, 2],
  accuracy(hwm5)[, 2]
))
rownames(RMSE) <- c(
    "Doanh thu trên 1 đồng chi phí sản xuất",
    "Doanh thu trên 1 đồng vốn kinh doanh",
    "Doanh lợi theo chi phí sản xuất",
    "Doanh lợi theo doanh thu",
    "Doanh lợi theo vốn kinh doanh"
  )
paged_table(RMSE)