LCD-TV Sales Data

The size profile of the sales data clearly indicates an exponential growth in sales during the period 2003–20 07. However, the plot of the ISRP profile of r reveals the true scenario. It is evident from the ISRP plot of r, that the ISRP is a decreasing function with respect to time, which indicates that exponential growth is not an appropriate model to choose. Rather it suggests that the exponential model with an intrinsic growth rate, which is a decreasing function of time, would be a better choice. So, for the analysis, we fit exponential, linear, and nonmonotonic functions for r(t). In particular, exponential decay of the form (r(t) = r0 exp[−ct]) has been found to be more appropriate than the linear (r(t) = r0 t) and nonmonotonic (r(t) = r0 exp[−ct] t^(a) ) choice of r(t) based on the AIC. AIC of exponential decay, linear, and nonmonotonic functions are -27.85451, -14.5925, and -26.9082, respectively. The model with the smallest AIC is considered as the best model. From a technical point of view, since the difference between the AIC values between r(t) = r0 exp[−ct] t^(a) and r(t) = r0 exp[−ct] (Gompertz model) is less than the threshold value, both models perform similarly. However, from the principle of parsimony, we selected the simpler model with a lesser number of parameters. Hence, the Gompertz model is the appropriate choice for the given data with the difference in AIC being greater than 10 in comparison to the linear model. We also fit all three models, namely, exponential, Gompertz, and logistic, and the respective AIC values for these models are 295.0443, 221.4758, and 247.4407, respectively. The smallest AIC for the Gompertz model indicates the preference of this model over the other two. The models are also compared using the root mean squared error, and the conclusion remains the same.

# The data LCD-TV: Trappey and Wu 2008

x = c( 26.900,  63.100, 164.100,  407.500, 787.500, 

       1194.500,  1603.900, 2178.600, 2993.600,  4059.600, 

       5432.200,  7101.652,  8883.652,  10896.652,  

       13379.652, 16097.652,  18563.652,  21794.652)

n = length(x)

t=0:(length(x)-1)

plot(t, x, lwd=2, col = "blue", cex=2, type = "p")  # plot of size profile

RGR = log(x[2:n]/x[1:(n-1)])     # plot of RGR profile

data = data.frame(t, x)

# Fitting of Gompertz Model

GompertzModel = function(t, x0, r0, alpha){

  x0*exp(r0*(1-exp(-alpha*t))/alpha)

}

par(mar = c(5.15, 5.15, 3.8, 1.8))

# Choosing the initial values of the parameters

initial_fit = nls(RGR ~ r0 * exp(-alpha*t[1:(n-1)]), start = list(r0=0.8, alpha=1))

r0_initial = coef(initial_fit)[1]; r0_initial

alpha_initial = coef(initial_fit)[2]; alpha_initial

plot(t[1:(n-1)], RGR, col = "blue", lwd=2, xlab="Time", 

     ylab = expression(bold(widehat(r[j](Delta~t)))), cex.lab = 1.5)

lines(t[1:(n-1)], fitted.values(initial_fit), col = "red", lwd=3, lty=2)

AIC(initial_fit)

abline(lm(RGR ~ t[1:(n-1)]), lwd = 3, col="magenta", lty=3)

AIC(lm(RGR ~ t[1:(n-1)]))

t1 = c(0.1, t[2:(n-1)])

new_fit = nls(RGR ~ r0 * (t1^c)*exp(-alpha*t1), 

              start = list(r0=1, alpha=1, c = 1))

AIC(new_fit)

lines(t1, fitted.values(new_fit), col = "green", lwd=3, lty=4)

legend("topright",

       legend = c(expression(paste(r=r[0]*e^{-ct})), expression(paste(r=r[0]*t)),

                  expression(paste(r=r[0]*e^{-ct}*t^{a}))), 

       lwd=c(3,3,3), col = c("red", "magenta", "green"), bty = "n",

       lty = c(2,3,4))

axis(1,font = 2); axis(2,font = 2)

# Fitting of the Gompertz model

gom_fit = nls(x ~ GompertzModel(t, x0, r0, alpha), 

              start = list(x0=x[1], r0 = r0_initial, alpha = alpha_initial))

est_x0 = coef(gom_fit)[1]

est_r0 = coef(gom_fit)[2]

est_alpha = coef(gom_fit)[3]

final_size = est_x0*exp(est_r0/est_alpha) # asymptotic size

print(final_size)

plot(t, x, lwd=2, cex=1.5, col = "blue", xlim = c(0,50), 

     ylim=c(0,final_size+20000), xlab = "Time", ylab="Size ",

     cex.lab=1.5, pch = 19)  # plot of size profile

curve((est_x0*exp(est_r0*(1-exp(-est_alpha*x))/est_alpha)), 

      lwd = 4, col = "green", add = TRUE, lty = 3)

# Prediction interval for the asymptotic size using delta method

library(car)

deltaMethod(gom_fit, "x0*exp(r0/alpha)")

AIC(gom_fit)

#x = as.numeric(readClipboard()) # load the complete data

t_vals = 0:(length(x)-1)

conf_band = matrix(data = NA, nrow = length(t_vals), ncol=2)

for (i in 1:length(t_vals)) {

  time = t_vals[i]

  m = deltaMethod(gom_fit, "x0*exp(r0*(1-exp(-alpha*time))/alpha)")

  conf_band[i,1] = m$`2.5 %`

  conf_band[i,2] = m$`97.5 %`

}

lines(t_vals, conf_band[,1], col = "red", lwd=3, lty = 4)

lines(t_vals, conf_band[,2], col = "red", lwd=3, lty = 4)

lines(t_vals, x, lwd=3, col = "grey", lty = 5)

# Fitting of exponential model

exp_fit = nls(x ~ x0*exp(r*t), start = list(x0 = x[1], r=0.2))

curve(coef(exp_fit)[1]*exp(coef(exp_fit)[2]*x), col = "magenta", 

      lwd=3, add = TRUE, lty = 4)

AIC(exp_fit)

# Fitting of logistic Model

LogisticModel = function(t, x0, r, K){

  K/(1+(K/x0 - 1)*exp(-r*t))

}

log_fit = nls(x ~ LogisticModel(t, x0, r, K), 

              start = list(x0 = x[1], r = 1, K = 10000) )

est_x0 = coef(log_fit)[1]

est_r = coef(log_fit)[2]

est_K = coef(log_fit)[3]

curve(est_K/(1+(est_K/est_x0 - 1)*exp(-est_r*x)), 

      lwd = 3, col = "red", add = TRUE, lty = 5)

AIC(log_fit)

legend("topleft", legend = c("Exponential", "Gompertz", "Logistic"),

       col = c("magenta", "green", "red"), lwd=c(3,3,3), bty="n",

       lty = c(4,3,5))

#points(t,x, col = "blue", lwd=2, cex=3)

print(c(AIC(exp_fit), AIC(gom_fit), AIC(log_fit)))

mean(residuals(exp_fit)^2)  # Exponential Fit

mean(residuals(gom_fit)^2)  # Gompertz Fit

mean(residuals(log_fit)^2)  # Logistic Fit

axis(1,font = 2); axis(2,font = 2)