Analytics for Observational Data - Alfie & Kiran
library("readr")
library("lattice")
library("ggbiplot")
## Warning: package 'ggbiplot' was built under R version 4.3.3
## Loading required package: ggplot2
data <- read_csv("C:/Users/Kiran/Desktop/AfOD Assessment Data - Monthly - CSV version.csv")
## Rows: 484 Columns: 57
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): Date
## dbl (56): Number of purchases ('000s) - Credit card, Value of purchases ($ m...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# Subset data
start_date <- "31/07/1998"
indexstart <- which(data==start_date)
hist_data <- data[indexstart:nrow(data),]
#subset hist data without date, and na.omit columns
hist_less_date <- hist_data[,c(2:57)]
hist_less_date_na <- hist_less_date[ , colSums(is.na(hist_less_date)) == 0]
# Quite different means and stdev so scale
colMeans(hist_less_date_na)
## Number of purchases ('000s) - Credit card
## 1.630043e+05
## Value of purchases ($ million) - Credit card
## 1.987448e+04
## Total number of transactions ('000s)
## 1.652394e+05
## Total value of transactions ($ million)
## 2.063643e+04
## Repayments ($ million)
## 2.096639e+04
## Credit; Owner-occupier housing; Seasonally adjusted ($ billion)
## 8.132821e+02
## Credit; Investor housing; Seasonally adjusted ($ billion)
## 4.108340e+02
## Credit; Other personal; Seasonally adjusted ($ billion)
## 1.278981e+02
## Short-term non-government securities issued in Australia – Total ($ billion)
## 2.745324e+02
## Long-term non-government securities issued in Australia – Total ($ billion)
## 3.943111e+02
## Consumer price index
## 9.961235e+01
## Real GDP ($ million)
## 4.898834e+05
## Household consumption ($ million)
## 2.466064e+05
## Business investment ($ million)
## 5.871531e+04
## Public demand ($ million)
## 1.151126e+05
## Changes in inventories ($ million)
## 8.586049e+02
## Gross national expenditure ($ million)
## 4.563720e+05
## Real household disposable income ($ million)
## 2.892292e+05
## Household saving ratio (per cent)
## 4.754012e+00
## Retail sales, seasonally adjusted ($ million)
## 2.200054e+04
## Private dwelling approvals (000's)
## 1.466636e+01
## Business conditions (Percentage)
## 3.675309e+00
## Year-ended wage growth (Percentage)
## 3.113580e+00
## Non-farm labour productivity per hour (Index)
## 9.278951e+01
## Labour force ('000s)
## 1.187602e+04
## Employment ('000s)
## 1.124471e+04
## Unemployment ('000s)
## 6.313120e+02
## Hours worked ('000s)
## 1.575964e+06
## Job vacancies ('000s, Quarterly)
## 1.927642e+02
## Exports ($ million)
## 1.311810e+05
## Imports ($ million)
## 9.895439e+04
## Trade balance as a per cent of output ($ million)
## -6.944444e-02
## Commodity prices – A$
## 6.007068e+01
## Rural commodity prices – A$ (Index, 2023/24=100)
## 7.996235e+01
## Non-rural commodity prices – A$ (Index, 2023/24=100)
## 5.802870e+01
## Base metals prices – A$ (Index, 2023/24=100)
## 6.694722e+01
## Bulk commodities prices – A$ (Index, 2023/24=100)
## 5.590463e+01
## Gross total foreign assets ($ billion)
## 1.782159e+03
## Gross total foreign liabilities ($ billion)
## 2.498210e+03
apply(hist_less_date_na, 2, sd)
## Number of purchases ('000s) - Credit card
## 8.175011e+04
## Value of purchases ($ million) - Credit card
## 8.871138e+03
## Total number of transactions ('000s)
## 8.107279e+04
## Total value of transactions ($ million)
## 8.751216e+03
## Repayments ($ million)
## 8.926263e+03
## Credit; Owner-occupier housing; Seasonally adjusted ($ billion)
## 4.202603e+02
## Credit; Investor housing; Seasonally adjusted ($ billion)
## 2.189344e+02
## Credit; Other personal; Seasonally adjusted ($ billion)
## 2.922552e+01
## Short-term non-government securities issued in Australia – Total ($ billion)
## 7.887597e+01
## Long-term non-government securities issued in Australia – Total ($ billion)
## 2.009569e+02
## Consumer price index
## 1.991911e+01
## Real GDP ($ million)
## 1.004040e+05
## Household consumption ($ million)
## 5.238077e+04
## Business investment ($ million)
## 1.768176e+04
## Public demand ($ million)
## 3.218609e+04
## Changes in inventories ($ million)
## 2.126572e+03
## Gross national expenditure ($ million)
## 1.047385e+05
## Real household disposable income ($ million)
## 6.995643e+04
## Household saving ratio (per cent)
## 4.634353e+00
## Retail sales, seasonally adjusted ($ million)
## 7.600891e+03
## Private dwelling approvals (000's)
## 2.708491e+00
## Business conditions (Percentage)
## 9.215565e+00
## Year-ended wage growth (Percentage)
## 7.891294e-01
## Non-farm labour productivity per hour (Index)
## 7.343586e+00
## Labour force ('000s)
## 1.689827e+03
## Employment ('000s)
## 1.666824e+03
## Unemployment ('000s)
## 9.844474e+01
## Hours worked ('000s)
## 1.992369e+05
## Job vacancies ('000s, Quarterly)
## 9.477911e+01
## Exports ($ million)
## 3.391096e+04
## Imports ($ million)
## 3.466334e+04
## Trade balance as a per cent of output ($ million)
## 2.321036e+00
## Commodity prices – A$
## 2.451211e+01
## Rural commodity prices – A$ (Index, 2023/24=100)
## 2.181103e+01
## Non-rural commodity prices – A$ (Index, 2023/24=100)
## 2.544977e+01
## Base metals prices – A$ (Index, 2023/24=100)
## 2.180971e+01
## Bulk commodities prices – A$ (Index, 2023/24=100)
## 2.857109e+01
## Gross total foreign assets ($ billion)
## 1.166146e+03
## Gross total foreign liabilities ($ billion)
## 1.345432e+03
# PCA n>p so prcomp over princomp
hist_pr <- prcomp(hist_less_date_na,scale=TRUE, center=TRUE)
summary(hist_pr)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 5.4163 1.7664 1.46940 1.28274 0.98504 0.73997 0.54154
## Proportion of Variance 0.7522 0.0800 0.05536 0.04219 0.02488 0.01404 0.00752
## Cumulative Proportion 0.7522 0.8322 0.88758 0.92977 0.95465 0.96869 0.97621
## PC8 PC9 PC10 PC11 PC12 PC13 PC14
## Standard deviation 0.46561 0.43433 0.33960 0.33419 0.30346 0.24638 0.1971
## Proportion of Variance 0.00556 0.00484 0.00296 0.00286 0.00236 0.00156 0.0010
## Cumulative Proportion 0.98177 0.98661 0.98957 0.99243 0.99479 0.99635 0.9973
## PC15 PC16 PC17 PC18 PC19 PC20 PC21
## Standard deviation 0.16422 0.14211 0.11757 0.10075 0.09906 0.08539 0.06082
## Proportion of Variance 0.00069 0.00052 0.00035 0.00026 0.00025 0.00019 0.00009
## Cumulative Proportion 0.99803 0.99855 0.99891 0.99917 0.99942 0.99961 0.99970
## PC22 PC23 PC24 PC25 PC26 PC27 PC28
## Standard deviation 0.05078 0.04523 0.03923 0.03432 0.03074 0.02905 0.02547
## Proportion of Variance 0.00007 0.00005 0.00004 0.00003 0.00002 0.00002 0.00002
## Cumulative Proportion 0.99977 0.99982 0.99986 0.99989 0.99991 0.99993 0.99995
## PC29 PC30 PC31 PC32 PC33 PC34
## Standard deviation 0.02408 0.02185 0.02038 0.01539 0.009917 0.007193
## Proportion of Variance 0.00001 0.00001 0.00001 0.00001 0.000000 0.000000
## Cumulative Proportion 0.99997 0.99998 0.99999 1.00000 1.000000 1.000000
## PC35 PC36 PC37 PC38 PC39
## Standard deviation 0.004732 0.003918 0.002021 0.0005288 1.907e-05
## Proportion of Variance 0.000000 0.000000 0.000000 0.0000000 0.000e+00
## Cumulative Proportion 1.000000 1.000000 1.000000 1.0000000 1.000e+00
hist_pr$loadings[,c(1:7)]
## NULL
# Screeplot
screeplot(hist_pr,type="lines",main="Screeplot of Variables 31/07/1998 to 30/06/2025")
abline(h=1,col="red")
# Test prediction against data (?) removed cuz huge output
# predict(hist_pr)
# dotplot PC1
load <- hist_pr$rotation
sorted.loadings <- load[order(load[, 1]), 1]
dotplot(sorted.loadings, main = "Loadings Plot for PC1", xlab = "Variable Loadings", col = "red")
sig.loadings<-load[order(abs(load[,1]),decreasing=TRUE),1]
head(sig.loadings,n=10)
## Real GDP ($ million)
## 0.1839603
## Credit; Owner-occupier housing; Seasonally adjusted ($ billion)
## 0.1837493
## Employment ('000s)
## 0.1836552
## Labour force ('000s)
## 0.1836527
## Gross national expenditure ($ million)
## 0.1835925
## Consumer price index
## 0.1829351
## Retail sales, seasonally adjusted ($ million)
## 0.1828393
## Household consumption ($ million)
## 0.1826734
## Credit; Investor housing; Seasonally adjusted ($ billion)
## 0.1825302
## Hours worked ('000s)
## 0.1824908
# dotplot PC2
sorted.loadings <- load[order(load[, 2]), 2]
dotplot(sorted.loadings, main = "Loadings Plot for PC2", xlab = "Variable Loadings", col = "red")
sig.loadings<-load[order(abs(load[,2]),decreasing=TRUE),2]
head(sig.loadings,n=10)
## Unemployment ('000s)
## 0.5168927
## Year-ended wage growth (Percentage)
## -0.4709652
## Private dwelling approvals (000's)
## 0.3490416
## Short-term non-government securities issued in Australia – Total ($ billion)
## -0.2303464
## Base metals prices – A$ (Index, 2023/24=100)
## -0.2263336
## Household saving ratio (per cent)
## 0.2250239
## Non-rural commodity prices – A$ (Index, 2023/24=100)
## -0.2054285
## Job vacancies ('000s, Quarterly)
## -0.1988841
## Commodity prices – A$
## -0.1984020
## Bulk commodities prices – A$ (Index, 2023/24=100)
## -0.1752351
# dotplot PC3
sorted.loadings <- load[order(load[, 3]), 3]
dotplot(sorted.loadings, main = "Loadings Plot for PC3", xlab = "Variable Loadings", col = "red")
sig.loadings<-load[order(abs(load[,3]),decreasing=TRUE),3]
head(sig.loadings,n=10)
## Short-term non-government securities issued in Australia – Total ($ billion)
## 0.5000829
## Changes in inventories ($ million)
## -0.4444720
## Business conditions (Percentage)
## -0.4138226
## Credit; Other personal; Seasonally adjusted ($ billion)
## 0.2835186
## Business investment ($ million)
## 0.2546395
## Job vacancies ('000s, Quarterly)
## -0.2097188
## Rural commodity prices – A$ (Index, 2023/24=100)
## -0.1725369
## Trade balance as a per cent of output ($ million)
## -0.1571910
## Private dwelling approvals (000's)
## -0.1415974
## Year-ended wage growth (Percentage)
## 0.1376904
# dotplot PC4
sorted.loadings <- load[order(load[, 4]), 4]
dotplot(sorted.loadings, main = "Loadings Plot for PC4", xlab = "Variable Loadings", col = "red")
sig.loadings<-load[order(abs(load[,4]),decreasing=TRUE),4]
head(sig.loadings,n=10)
## Business conditions (Percentage)
## -0.4819552
## Trade balance as a per cent of output ($ million)
## 0.4268618
## Household saving ratio (per cent)
## 0.3694605
## Private dwelling approvals (000's)
## -0.3410940
## Bulk commodities prices – A$ (Index, 2023/24=100)
## 0.2449054
## Short-term non-government securities issued in Australia – Total ($ billion)
## -0.1941579
## Credit; Other personal; Seasonally adjusted ($ billion)
## -0.1890303
## Commodity prices – A$
## 0.1812075
## Non-rural commodity prices – A$ (Index, 2023/24=100)
## 0.1783561
## Rural commodity prices – A$ (Index, 2023/24=100)
## 0.1619534
g <- ggbiplot(hist_pr, choices=c(1,2), obs.scale = 1, var.scale = 1, ellipse = TRUE, circle = TRUE)
g <- g + scale_color_discrete(name = '')
g <- g + theme(legend.direction = 'horizontal',
legend.position = 'top')
print(g+ylim(-10,10))
# Correlated
image(cor(hist_less_date_na))
library("fable")
## Warning: package 'fable' was built under R version 4.3.3
## Loading required package: fabletools
## Registered S3 method overwritten by 'tsibble':
## method from
## as_tibble.grouped_df dplyr
library("tsibble")
## Warning: package 'tsibble' was built under R version 4.3.3
##
## Attaching package: 'tsibble'
## The following objects are masked from 'package:base':
##
## intersect, setdiff, union
library("feasts")
library("forecast")
## Warning: package 'forecast' was built under R version 4.3.3
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
CPI_date<-hist_data[,c(1,28)]
plot(CPI_date$`Consumer price index`,xlab="Month", ylab="CPI")
# Does not look stationary, take first difference
CPI_ts <- ts(CPI_date[,2],start=1,frequency=1)
ggtsdisplay(diff(CPI_ts,main=""))
CPI_tsibble <- as_tsibble(CPI_ts)
# Fit multiple models
models <- model(CPI_tsibble,
ETS = ETS(value),
ARIMA = ARIMA(value)
)
# Generate forecasts
fc <- forecast(models, h = 36)
# Plot past + forecasts together
autoplot(fc, CPI_tsibble, level = NULL) +
ggplot2::labs(title = "CPI Forecasts: ETS vs ARIMA",
y = "CPI", x = "Date") +
ggplot2::theme_minimal()
report(models[1])
## Series: value
## Model: ETS(A,Ad,N)
## Smoothing parameters:
## alpha = 0.9759875
## beta = 0.9749844
## phi = 0.8798708
##
## Initial states:
## l[0] b[0]
## 67.40729 -0.1391341
##
## sigma^2: 0.0231
##
## AIC AICc BIC
## 658.9458 659.2108 681.6303
report(models[2])
## Series: value
## Model: ARIMA(3,2,3)
##
## Coefficients:
## ar1 ar2 ar3 ma1 ma2 ma3
## -0.0290 -0.0683 0.2083 -0.1648 0.0259 -0.7482
## s.e. 0.1038 0.0879 0.0979 0.0799 0.0659 0.0738
##
## sigma^2 estimated as 0.01865: log likelihood=186.53
## AIC=-359.06 AICc=-358.7 BIC=-332.64
res1 <- augment(models[1])
autoplot(res1,.resid,main="Residuals over time")
## Warning in geom_line(...): Ignoring unknown parameters: `main`
features(res1,.resid, ljung_box, lag = 20)
## # A tibble: 1 × 3
## .model lb_stat lb_pvalue
## <chr> <dbl> <dbl>
## 1 ETS 44.5 0.00131
res2 <- augment(models[2])
autoplot(res2,.resid,main="Residuals over time")
## Warning in geom_line(...): Ignoring unknown parameters: `main`
features(res2,.resid, ljung_box, lag = 20)
## # A tibble: 1 × 3
## .model lb_stat lb_pvalue
## <chr> <dbl> <dbl>
## 1 ARIMA 4.24 1.00