• What Is Data Cleaning?

• Data Wrangling (aka Data Munging)

• Outlier Detection

  • Statistical and likelihood approaches
  • Depth-based, deviation-based, distance-based, density-based approaches
  • High-dimensional approaches

• Imputation of Missing Values

• Data Cleaning of Financial Time Series

  • Robust model fitting
  • Outlier detection
  • Imputation of missing values

Consider a simple linear regression:

set.seed(42)
x <- runif(40)
y <- 1 + 3*x + 0.1*rnorm(20)
res <- lm(y ~ x)
res

R>> 
R>> Call:
R>> lm(formula = y ~ x)
R>> 
R>> Coefficients:
R>> (Intercept)            x  
R>>      0.9818       2.9849

# plot
plot(x, y, main = "Regression without outliers")
abline(res, col = "red")

Now let’s see what happens with a couple of outliers:

# add two outliers
y[35] <- y[35] + 3
y[23] <- y[23] - 4
res <- lm(y ~ x)
res

R>> 
R>> Call:
R>> lm(formula = y ~ x)
R>> 
R>> Coefficients:
R>> (Intercept)            x  
R>>       1.506        2.054

We can see that the intercept is estimated as 1.5 instead of 1 and the slope is estimated as 2 instead of 3!

The plot shows how detrimental the effect is:

# plot
plot(x, y, main = "Regression with outliers")
abline(res, col = "red")

Let’s first create simple data with one large outlier:

# create data with one large outlier
data_with_oulier <- c(1:10, 300)

# mean of clean data
mean(data_with_oulier[1:10])

R>> [1] 5.5

# mean of data with outlier
mean(data_with_oulier)

R>> [1] 32.27273

Trimming:

The mean() function in R has a trim argument so that you can easily get trimmed means. To be exaxt, trim denotes the fraction (0 to 0.5) of observations to be trimmed from each end of the data before the mean is computed.

# trimming 5% of the data is not enough to get rid of the outlier
mean(data_with_oulier, trim = .05)

R>> [1] 32.27273

# trimming 10% does the job
mean(data_with_oulier, trim = .10)

R>> [1] 6

# trimming removes both extremes, so it's like doing:
mean(data_with_oulier[2:10])

R>> [1] 6

Winsorization:

We will now define a function to Winsorize:

winsorize <- function (x, fraction = .05) {
  lim <- quantile(x, probs = c(fraction, 1-fraction))
  x[x < lim[1]] <- lim[1]
  x[x > lim[2]] <- lim[2]
  return(x)
}

winsorize(data_with_oulier, fraction = .10)

R>>  [1]  2  2  3  4  5  6  7  8  9 10 10

mean(winsorize(data_with_oulier, fraction = .10))

R>> [1] 6

For this study on detecting outliers in R, we will follow the example in RPubs. We will be using the csv file “PregnancyDuration.csv” available for download at www.wiley.com/go/datasmart.

pregnancy.df <- read.csv("PregnancyDuration.csv")
str(pregnancy.df)

R>> 'data.frame':   1000 obs. of  1 variable:
R>>  $ GestationDays: int  349 278 266 265 269 263 278 257 268 260 ...

head(pregnancy.df)

R>>   GestationDays
R>> 1           349
R>> 2           278
R>> 3           266
R>> 4           265
R>> 5           269
R>> 6           263

summary(pregnancy.df)

R>>  GestationDays  
R>>  Min.   :240.0  
R>>  1st Qu.:260.0  
R>>  Median :267.0  
R>>  Mean   :266.6  
R>>  3rd Qu.:272.0  
R>>  Max.   :349.0

From the output, we see the minimum gestation days is 240 days or 8 months, while the maximum gestation days is 349 days or about 11 months. Therefore, the range of days between gestation is 109 days or about 3 months. This may suggest lots of variability in gestation days. Plus, the 25th percentile of data is 260 days or about 8,66 months and the 75th percentile is 272 or about 9,06 months. Finally, the median gestation days is 267 and the “typical” average gestation days is about 267 days or exactly 9 months. Let’s check the histogram of gestation days.

Let’s plot the histogram to visualize the data better:

library(ggplot2)

# Histogram of gestation days
qplot(pregnancy.df$GestationDays, geom = "histogram", fill = I("cyan"), col = I("black"),
      main = "Histogram of Gestation Days", xlab = "Gestation Days")

The output has a peaked bar with gestation days at about 272 days. Most of the values are bunched up on the left side of the histogram, with few values spreading along the right tail. This tells us that the data points are not normally distributed and it has an asymmetric distribution (note that the mean and the median occur at different points). In addition, the plot reveals the presence of an outlier in the lower right bound of the plot. This outlier is completely apart from other observations.

Using the histogram is a good way to assess the shape and spread of the data and to identify any potential outliers. Similarly, the boxplot is also vital to evaluate the spread of the data. Let’s go back to the summary and get a good understanding of skewness.

Interquartile range: In the summary results, the interquartile range is equal to 272 minus 260, IQR = 3rd qu - 1st qu. Alternatively, you can call the built-in IQR() function on the GestationDays column to calculate the IQR:

summary(pregnancy.df)

R>>  GestationDays  
R>>  Min.   :240.0  
R>>  1st Qu.:260.0  
R>>  Median :267.0  
R>>  Mean   :266.6  
R>>  3rd Qu.:272.0  
R>>  Max.   :349.0

pregnancy_IQR <- IQR(pregnancy.df$GestationDays)
pregnancy_IQR

R>> [1] 12

# we can also find the quartiles directly
quantile(pregnancy.df$GestationDays, c(0.25, 0.75))

R>> 25% 75% 
R>> 260 272

We can then compute the lower and upper Tukey fences, thanks to John Wilder Tukey. John Wilder Tukey was an American mathematician best known for the development of the FFT algorithm and boxplot. The Tukey range test, the Tukey lambda distribution, the Tukey test of additivity, and the Teichmuller-Tukey lemma all bear his name.

lower_inner_fence <- 260 - 1.5 * pregnancy_IQR
lower_inner_fence

R>> [1] 242

upper_inner_fence <- 272 + 1.5 * pregnancy_IQR
upper_inner_fence

R>> [1] 290

Let’s see the points that violate the Tukey fences:

# values above the upper Tukey fence
pregnancy.df$GestationDays[pregnancy.df$GestationDays > upper_inner_fence]

R>> [1] 349 292 295 291 297 303 293 296

# values below the lower Tukey fence
pregnancy.df$GestationDays[pregnancy.df$GestationDays < lower_inner_fence]

R>> [1] 241 240

Two points are below the lower Tukey fence and eight values are above the upper Tukey fence.

Boxplot:

We are now ready to plot a boxplot using the function boxplot(). The boxplot() function will graph the median, first quartile and third quartile, Tukey fences, and any outliers.

boxplot(pregnancy.df$GestationDays)

The tukey fences can be modified to be “outlier” fences by changing the range flag in the boxplot call(it defaults to 1.5 times the IQR). If you set range = 3, then the tukey fences are drawn at the last point inside three times the IQR instead.

boxplot(pregnancy.df$GestationDays, range = 3)

We can improve the boxplot by using ggplot2. You can draw a boxplot and show the colored outliers.

# boxplot using ggplot2 with colored outliers
ggplot(data = pregnancy.df, mapping = aes(x = "", y = GestationDays)) +
  geom_boxplot(outlier.color = "red", outlier.shape = 1) +
  ggtitle("Boxplot of The Pregnancy Duration Data", subtitle = "With red outliers in Tukey fences using 1.5 times the IQR") +
  xlab("") +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(plot.subtitle = element_text(hjust = 0.5, face = "italic"))

In the output result, the tiny red cirles are outliers. You have 8 outliers above the box and 2 outliers below the box. Let’s expand our plot by adding the colored mean to the boxplot.

We will explore the detection of outliers in a multivariate setting based on Grubbs’ test using Mahalanobis distance.

Let’s generate multivariate data and add some outliers:

library(mvtnorm)  # for multivariate Gaussian and Student t generation
library(ggplot2)

# create a mean vector and covariance matrix
p <- 2   # dimension: we use 2 so that we can plot and visualize
n <- 30  # number of points
set.seed(42)
mu <- rep(0, p)
U <- t(rmvnorm(n = ceiling(0.1*p), sigma = 0.1*diag(p)))
Sigma <- U %*% t(U) + 0.1*diag(p)

# generate n x p Gaussian data matrix and outliers
X <- rmvnorm(n = n, mean = mu, sigma = Sigma)
# add four outliers (in the first four data points)
X[1, 2] <- 8
X[2, 2] <- -9
X[3, 2] <- 7
X[4, 2] <- -7
X_clean    <- X[-c(1:4), ]
X_outliers <- X[ c(1:4), ]

We can now plot the generated points together with the ellipsoid of the clean data (black) and the ellipsoid including the outliers (red):

# 2D plots from the first two components of data matrix
ggplot(data = as.data.frame(X_clean)) +
  geom_point(aes(V1, V2)) +
  stat_ellipse(aes(V1, V2)) +
  geom_point(data = as.data.frame(X_outliers), aes(V1, V2), col = "red") +
  stat_ellipse(data = as.data.frame(X), aes(V1, V2), col = "red") +
  ggtitle("Gaussian data with four outliers") + xlab(NULL) + ylab(NULL)

Indeed, the estimated covariance matrix with outliers (red ellipsoid) is severely affected by the outliers.

Let’s try to detect them based on the Mahalanobis distance:

Sigma <- cov(X)
inv_Sigma <- solve(Sigma)
Mahalanobis_distance <- colSums(t(X) * (inv_Sigma %*% t(X)))
qplot(1:n, Mahalanobis_distance, main = "Mahalanobis distances of datapoints", xlab = "", ylab = "")

We can see that the four outliers are clearly identified, but other points come close too. To be more accurate, after we remove the outliers in the estimation of the covariance matrix, we should re-compute again the Mahalanobis distances:

# re-compute Sigma without the identified outliers
Sigma <- cov(X[-c(1:4), ])
inv_Sigma <- solve(Sigma)
Mahalanobis_distance <- colSums(t(X) * (inv_Sigma %*% t(X)))
qplot(1:n, Mahalanobis_distance, main = "Mahalanobis distances of datapoints", xlab = "", ylab = "")

Now the outliers are clearly identified. In practice, one can remove the outliers one by one.

In this session, we will illustrate the detrimental effects of outliers in the estimation of the mean vector \(\boldsymbol{\mu}\) and covariance matrix \(\boldsymbol{\Sigma}\).

The naive benchmark estimators are the sample mean and sample covariance matrix: \[\hat{\boldsymbol{\mu}} = \frac{1}{T}\sum_{t=1}^T \mathbf{x}_t\] and \[\hat{\boldsymbol{\Sigma}} = \frac{1}{T}\sum_{t=1}^T (\mathbf{x}_t - \hat{\boldsymbol{\mu}})(\mathbf{x}_t - \hat{\boldsymbol{\mu}})^T.\]

For the robust estimators, we will use the R package fitHeavyTail that assumes the data follows a heavy-tailed distribution (which is a way to model naturally outliers).

library(fitHeavyTail)  # for robust estimation
library(mvtnorm)  # for multivariate Gaussian and Student t generation
library(ggplot2)
library(reshape2)  # for melt()

# create a mean vector and covariance matrix
N <- 50  # dimension
set.seed(42)
mu <- runif(N)
U <- t(rmvnorm(n = ceiling(0.1*N), sigma = 0.1*diag(N)))
Sigma <- U %*% t(U) + 0.1*diag(N)
nu <- 4  # for heavy-tailed distribution
Sigma_scale <- (nu-2)/nu*Sigma   # Sigma = nu/(nu-2)*Sigma_scale

# loop over different values of the number of observations T
T_sweep <- round(seq(from = N+1, to = 5*N, length.out = 50))
MSE_mu_robust_T <- MSE_Sigma_robust_T <- MSE_sm_T <- MSE_scm_T <- NULL
for(T in T_sweep) {
  # generate data matrix
  #X <- rmvnorm(n = T, mean = mu, sigma = Sigma)  # Gaussian data
  X <- rmvt(n = T, delta = mu, sigma = Sigma_scale, df = nu)  # heavy-tailed data

  # naive estimation (sample estimates)
  mu_sm <- colMeans(X)
  Sigma_scm <- cov(X)

  # robust estimation against outliers and heavy tails
  moments_Studentt <- fit_mvt(X)
  mu_robust    <- moments_Studentt$mu
  Sigma_robust <- moments_Studentt$cov
  
  # compute MSEs
  MSE_sm_T  <- c(MSE_sm_T,  sum((mu_sm - mu)^2))
  MSE_scm_T <- c(MSE_scm_T, norm(Sigma_scm - Sigma, "F")^2)
  MSE_mu_robust_T    <- c(MSE_mu_robust_T,    sum((mu_robust - mu)^2))
  MSE_Sigma_robust_T <- c(MSE_Sigma_robust_T, norm(Sigma_robust - Sigma, "F")^2)
}

# MSE plots
MSE_Sigma_T <- cbind("SCM"     = MSE_scm_T, 
                     "robust"  = MSE_Sigma_robust_T)
MSE_mu_T <- cbind("sample mean" = MSE_sm_T, 
                  "robust"      = MSE_mu_robust_T)
rownames(MSE_mu_T) <- rownames(MSE_Sigma_T) <- T_sweep

ggplot(melt(MSE_mu_T), aes(x = Var1, y = value, col = Var2, shape = Var2)) + 
  geom_line() + geom_point() +
  theme(legend.title = element_blank()) +
  ggtitle(sprintf("MSE of the estimation of mu (N = %d)", N)) + xlab("T") + ylab("MSE")

ggplot(melt(MSE_Sigma_T), aes(x = Var1, y = value, col = Var2, shape = Var2)) + 
  geom_line() + geom_point() +
  theme(legend.title = element_blank()) +
  ggtitle(sprintf("MSE of the estimation of Sigma (N = %d)", N)) + xlab("T") + ylab("MSE")

In this session, we will use the package imputeTS to illustrate some simple univariate imputation methods:

library(imputeTS)

For this we first need to create an example input series with missing data:

# Create a short example time series with missing values
x <- ts(c(1, 2, 3, 4, 5, 6, 7, 8, NA, NA, 11, 12))

On this time series we can apply different imputation algorithms. We start with using na_mean(), which substitutes the NAs with mean values:

# Impute the missing values with na_mean()
na_mean(x)

R>> Time Series:
R>> Start = 1 
R>> End = 12 
R>> Frequency = 1 
R>>  [1]  1.0  2.0  3.0  4.0  5.0  6.0  7.0  8.0  5.9  5.9 11.0 12.0

Most of the functions also have additional options that provide further algorithms (of the same algorithm category). In the example below it can be seen that na_mean() can also be called with option = "median", which substitutes the NAs with median values:

# Impute the missing values with na_mean() using option median
na_mean(x, option = "median")

R>> Time Series:
R>> Start = 1 
R>> End = 12 
R>> Frequency = 1 
R>>  [1]  1.0  2.0  3.0  4.0  5.0  6.0  7.0  8.0  5.5  5.5 11.0 12.0

While na_interpolation() and all other imputation functions are used the same way, the results produced may be different. As can be seen below, for this series linear interpolation gives more reasonable results:

# Impute the missing values with na_interpolation()
na_interpolation(x)

R>> Time Series:
R>> Start = 1 
R>> End = 12 
R>> Frequency = 1 
R>>  [1]  1  2  3  4  5  6  7  8  9 10 11 12

For longer and more complex time series (with trend and seasonality) than in this example it is always a good idea to try na_kalman() and na.seadec(), since these functions very often produce the best results. These functions are called the same easy way as all other imputation functions. Consider the time series dataset tsAirgap provided in the package. We can use the functions plotNA.distribution() and plotNA.distributionBar() to visualize the distribution of missing values within a time series:

# Visualize the missing values in this time series
ggplot_na_distribution(tsAirgap)

We can also use the functino statsNA() to print summary stats about the distribution of missing values in a univariate time series:

statsNA(tsAirgap)

R>> [1] "Length of time series:"
R>> [1] 144
R>> [1] "-------------------------"
R>> [1] "Number of Missing Values:"
R>> [1] 13
R>> [1] "-------------------------"
R>> [1] "Percentage of Missing Values:"
R>> [1] "9.03%"
R>> [1] "-------------------------"
R>> [1] "Number of Gaps:"
R>> [1] 11
R>> [1] "-------------------------"
R>> [1] "Average Gap Size:"
R>> [1] 1.181818
R>> [1] "-------------------------"
R>> [1] "Stats for Bins"
R>> [1] "  Bin 1 (36 values from 1 to 36) :      4 NAs (11.1%)"
R>> [1] "  Bin 2 (36 values from 37 to 72) :      1 NAs (2.78%)"
R>> [1] "  Bin 3 (36 values from 73 to 108) :      5 NAs (13.9%)"
R>> [1] "  Bin 4 (36 values from 109 to 144) :      3 NAs (8.33%)"
R>> [1] "-------------------------"
R>> [1] "Longest NA gap (series of consecutive NAs)"
R>> [1] "3 in a row"
R>> [1] "-------------------------"
R>> [1] "Most frequent gap size (series of consecutive NA series)"
R>> [1] "1 NA in a row (occuring 10 times)"
R>> [1] "-------------------------"
R>> [1] "Gap size accounting for most NAs"
R>> [1] "1 NA in a row (occuring 10 times, making up for overall 10 NAs)"
R>> [1] "-------------------------"
R>> [1] "Overview NA series"
R>> [1] "  1 NA in a row: 10 times"
R>> [1] "  3 NA in a row: 1 times"

Here is a usage example for the na_kalman() function:

# Impute the missing values with na_kalman
x_imp <- na_kalman(tsAirgap)

# Code for visualization
ggplot_na_imputations(tsAirgap, x_imp, tsAirgapComplete)

As can be seen, na_kalman() provides really good results for this series, which contains a strong seasonality and a strong trend.

However, for financial time series, the package imputeTS does not seem very appropriate, especially for consecutive NAs. Let’s download some financial data and artificially generate NAs:

library(quantmod)
# download data from YahooFinance
y <- log(Ad(getSymbols("^GSPC", from = "2012-01-01", to = "2015-07-08", auto.assign = FALSE)))
T <- nrow(y)

# create NA's
miss_pct <- 0.2
n_miss <- floor(miss_pct*T)
idx_miss <- round(T/2) + 1:n_miss
y[idx_miss] <- NA

# plot
plot(y, main = "Log-prices of SP500 with a block of NAs")

Let’s impute with package imputeTS:

y_locf   <- na_locf(y)
y_linear <- na_interpolation(y, option = "linear")
y_spline <- na_interpolation(y, option = "spline")
y_stine  <- na_interpolation(y, option = "stine")
y_kalman <- na_kalman(y)

# plots
idx_miss_bis <- (min(idx_miss)-1):(max(idx_miss)+1)

{ plot(y, main = "Imputation - LOCF")
  lines(y_locf[idx_miss_bis], col = "blue", lwd = 2) }

{ plot(y, main = "Imputation - linear")
  lines(y_linear[idx_miss_bis], col = "blue", lwd = 2) }

{ plot(y, main = "Imputation - spline")
  lines(y_spline[idx_miss_bis], col = "blue", lwd = 2) }

{ plot(y, main = "Imputation - stine")
  lines(y_stine[idx_miss_bis], col = "blue", lwd = 2) }

{ plot(y, main = "Imputation - Kalman")
  lines(y_kalman[idx_miss_bis], col = "blue", lwd = 2) }

For time series, it is better to use the package imputeFin:

library(imputeFin)

y_ar1 <- impute_AR1_Gaussian(y)

# plot
{ plot(y, main = "Imputation - AR(1)")
  lines(y_ar1[idx_miss_bis], col = "blue", lwd = 2) }

In this R session we will follow the complete procedure of data cleaning of a financial time series, namely, the S&P 500. We will first identify and remove the outliers and then impute the missing values (and removed outliers). We will use the R package imputeFin.

Let’s start by loading some market data and adding some outliers and missing values:

library(xts)
library(quantmod)

# download data from YahooFinance
y_orig <- log(Ad(getSymbols("^GSPC", from = "2012-01-01", to = "2015-07-08", auto.assign = FALSE)))
y <- y_orig
T <- nrow(y)

# create some outliers
idx_outliers <- round(c(T/4, 1.5*T/4, 3*T/4))
y[idx_outliers] <- 1.05*y[idx_outliers]

# create NA's
miss_pct <- 0.2
n_miss <- floor(miss_pct*T)
idx_miss <- round(T/2) + 1:n_miss
y[idx_miss] <- NA

# plot
plot(y, main = "Log-prices of SP500 with outliers and NAs")

The first step is to identify the outliers. We will fit the time series to some parametric model and then identify the samples that deviate from the model significantly. In particular, we will consider the AR(1) model: \[y_{t}=\phi_0+\phi_1 y_{t-1}+\varepsilon_{t},\] where the innovation \(\varepsilon_{t}\) follows a Student’s \(t\)-distribution \(\varepsilon_t\sim t\left(0,\sigma^{2},\nu\right),\) and then we will identify the innovation samples that deviate more than \(3\sigma-5\sigma\) (so basically Grubb’s test on the innovations of the parametric model).

It’s instructive to see the effect of the outliers if we don’t use a robust estimation method (i.e., assuming a Gaussian distribution for the innovation):

library(imputeFin)  # package to fit models with NAs and robust to outliers

# if we use a Gaussian (nonrobust)
res_y_orig_Gaussian <- fit_AR1_Gaussian(y_orig)

R>> GSPC.Adjusted: 0 inner missing values and 0 outliers detected.

res_y_Gaussian      <- fit_AR1_Gaussian(y, zero_mean = FALSE)

R>> GSPC.Adjusted: 176 inner missing values and 0 outliers detected.

cbind(res_y_orig_Gaussian, res_y_Gaussian)

R>>                res_y_orig_Gaussian res_y_Gaussian
R>> phi0           0.01486504          0.1612201     
R>> phi1           0.9980751           0.9783727     
R>> sigma2         5.489247e-05        0.001205647   
R>> index_miss     NULL                Integer,176   
R>> index_outliers NULL                NULL

The estimation of the parameters \(\sigma^2\) and \(\phi_0\) have been severely affected (\(\phi_1\) not so much).

Let’s try a robust estimation (i.e., assuming a Student’s \(t\)-distribution for the innovation):

res_y_orig_Studentt <- fit_AR1_t(y_orig, zero_mean = FALSE)

R>> GSPC.Adjusted: 0 inner missing values and 0 outliers detected.

res_y_Studentt      <- fit_AR1_t(y, zero_mean = FALSE)

R>> GSPC.Adjusted: 176 inner missing values and 0 outliers detected.

cbind(res_y_orig_Studentt, res_y_Studentt)

R>>                res_y_orig_Studentt res_y_Studentt
R>> phi0           0.01090953          0.01049288    
R>> phi1           0.9986277           0.9986609     
R>> sigma2         3.703216e-05        2.741813e-05  
R>> nu             5.778218            2.2622        
R>> index_miss     NULL                Integer,176   
R>> index_outliers NULL                NULL

Much better!

Let’s now identify the outliers:

res_y_Studentt      <- fit_AR1_t(y, zero_mean = FALSE, remove_outliers = TRUE)

R>> GSPC.Adjusted: 179 inner missing values and 3 outliers detected.

cbind(res_y_orig_Studentt, res_y_Studentt)

R>>                res_y_orig_Studentt res_y_Studentt
R>> phi0           0.01090953          0.01307457    
R>> phi1           0.9986277           0.9983113     
R>> sigma2         3.703216e-05        3.937763e-05  
R>> nu             5.778218            5.950581      
R>> index_miss     NULL                Integer,179   
R>> index_outliers NULL                Integer,3

At this point, we are ready to do the imputation of the missing values:

y_imputed <- impute_AR1_t(y, remove_outliers = TRUE)

R>> GSPC.Adjusted: 179 missing values imputed and 3 outliers detected and corrected.

#plot(y_imputed, main = "Log-prices of SP500 after removing outliers and imputing NAs")

To summarize, let’s compare the original dirty time series with the clean one:

plot_imputed(y, title = "Log-prices of SP500 with outliers and NAs")

plot_imputed(y_imputed, title = "Log-prices of SP500 after removing outliers and imputing NAs")