An Alternative Data-Driven Prediction Approach Based on Real Option Theories

TL;DR

This paper introduces a novel prediction model for volatile time series data by integrating a time-varying Geometric Brownian Motion with real option pricing, outperforming traditional linear models in accuracy.

Contribution

The paper develops a flexible, non-linear prediction model based on real option theory, addressing volatility and non-linearity in time series forecasting.

Findings

Demonstrates improved prediction accuracy over traditional models.

Applicable across manufacturing, finance, and environmental datasets.

Shows competitive performance in real-world scenarios.

Abstract

This paper presents a new prediction model for time series data by integrating a time-varying Geometric Brownian Motion model with a pricing mechanism used in financial engineering. Typical time series models such as Auto-Regressive Integrated Moving Average assumes a linear correlation structure in time series data. When a stochastic process is highly volatile, such an assumption can be easily violated, leading to inaccurate predictions. We develop a new prediction model that can flexibly characterize a time-varying volatile process without assuming linearity. We formulate the prediction problem as an optimization problem with unequal overestimation and underestimation costs. Based on real option theories developed in finance, we solve the optimization problem and obtain a predicted value, which can minimize the expected prediction cost. We evaluate the proposed approach using multiple datasets obtained from real-life applications including manufacturing, finance, and environment. The numerical results demonstrate that the proposed model shows competitive prediction capability, compared with alternative approaches.

Tables

Table 1. TABLE I: CM Prediction Results for ten types of bumper beams with ω = 1 / 1.15 𝜔 1 1.15 \omega=1/1.15 in the Testing Set (The values in bold indicate the lowest prediction error for each product)

Weighted Mean Absolute Error (WMAE)
Product No.	ARIMA	Option Prediction	ARIMA-GARCH	CM Prediction
1	1196.58	411.97	2105.25	1161.56
2	332.36	127.90	168.72	151.24
3	119.35	105.49	107.43	194.69
4	1476.09	574.52	2185.30	936.94
5	1330.40	1299.00	1327.08	1797.74
6	542.33	64.38	63.24	42.90
7	357.38	24.14	497.54	92.17
8	1776.17	1339.11	3103.77	2520.75
9	1496.62	1305.49	2887.48	2475.71
10	1125.52	516.06	3278.77	1928.58

Table 2. TABLE II: CM Prediction Results for ten types of bumper beams with ω = 1 𝜔 1 \omega=1 in the Testing Set (The values in bold indicate the lowest prediction error for each product)

Weighted Mean Absolute Error (WMAE)
Product No.	ARIMA	Option Prediction	ARIMA-GARCH	CM Prediction
1	1193.29	537.94	2103.43	1335.79
2	328.39	134.67	164.49	168.93
3	111.78	118.39	103.40	223.89
4	1307.50	609.06	1920.83	1073.16
5	1178.90	1403.17	1200.81	2060.95
6	472.86	74.52	57.08	45.12
7	315.42	34.61	434.33	93.51
8	1691.28	1465.94	2737.90	2590.93
9	1453.32	1397.14	2564.50	2527.94
10	1008.59	569.58	2865.30	1932.47

Table 3. TABLE III: CM Prediction Results for ten types of bumper beams with ω = 1.15 𝜔 1.15 \omega=1.15 in the Testing Set (The values in bold indicate the lowest prediction error for each product)

Weighted Mean Absolute Error (WMAE)
Product No.	ARIMA	Option Prediction	ARIMA-GARCH	CM Prediction
1	1196.58	411.97	2105.25	1161.56
2	332.36	127.90	168.72	151.24
3	119.35	105.49	107.43	194.69
4	1476.09	574.52	2185.30	936.94
5	1330.40	1299.00	1327.08	1797.74
6	542.33	64.38	63.24	42.90
7	357.38	24.14	497.54	92.17
8	1776.17	1339.11	3103.77	2520.75
9	1496.62	1305.49	2887.48	2475.71
10	1125.52	516.06	3278.77	1928.58

Table 4. TABLE IV: Dow Jones Index Price Prediction Results in the Testing Set (The values in bold indicate the lowest prediction error for each testing period and weight)

Testing Period	Weight ($\omega$)	Method	WMAE	WMAPE
		ARIMA-GARCH	595.54	0.04998
	1/1.15	Option Prediction	50.40	0.0043
		ARIMA	628.69	0.0615
		ARIMA-GARCH	594.58	0.0499
Oct 2010- Mar 2011	1	Option Prediction	54.30	0.0046
		ARIMA	553.05	0.0541
		ARIMA-GARCH	682.72	0.0572
	1.15	Option Prediction	58.78	0.0050
		ARIMA	559.59	0.0548
		ARIMA-GARCH	360.39	0.0237
	1/1.15	Option Prediction	70.94	0.0046
		ARIMA	455.35	0.0309
		ARIMA-GARCH	317.34	0.0209
Aug 2013 - Dec 2013	1	Option Prediction	75.78	0.0049
		ARIMA	406.97	0.0276
		ARIMA-GARCH	321.71	0.0211
	1.15	Option Prediction	81.09	0.0052
		ARIMA	418.90	0.0284
		ARIMA-GARCH	338.71	0.0195
	1/1.15	Option Prediction	97.78	0.0056
		ARIMA	182.60	0.0109
		ARIMA-GARCH	319.72	0.0184
Oct 2014 - Mar 2015	1	Option Prediction	104.70	0.0060
		ARIMA	166.65	0.0099
		ARIMA-GARCH	348.25	0.0200
	1.15	Option Prediction	113.11	0.0065
		ARIMA	175.25	0.0104

Table 5. TABLE V: Wind Speed Prediction Results in the Testing Set (The values in bold indicate the lowest prediction error for each weight)

Weight ($\omega$)	Method	WMAE	WMAPE
	ARIMA-GARCH	3.364	0.402
1/1.15	Option Prediction	0.318	0.040
	ARIMA	8.73	0.957
	ARIMA-GARCH	3.31	0.380
1	Option Prediction	0.342	0.043
	ARIMA	8.73	0.96
	ARIMA-GARCH	3.74	0.422
1.15	Option Prediction	0.365	0.046
	ARIMA	10.03	1.100

Equations

$$dS(t)=\mu (t)S(t)dt+\sigma (t)S(t)dW(t),$$

$$S(T)=S(t)\>exp\Big{(}\int_{t}^{T}\Big{(}\mu(s)-\frac{1}{2}\sigma^{2}(s)\Big{)}ds+\int_{t}^{T}\sigma(s)dW(s)\Big{)},$$

$$\mathbb{E}(S(T)|\mathcal{F}_{t})=S(t)\>exp\big{(}\int_{t}^{T}\mu(s)ds\big{)}.$$

$$\displaystyle\displaystyle\min_{K\in R^{+}}\mathbb{E}\Big{[}P_{o}\max\{K-S(T),0\}+P_{u}\max\{S(T)-K,0\}|\mathcal{F}_{t}\Big{]}.$$

$$\max \{K-S(T),0\}=K-S(T)+\max \{S(T)-K,0\}.$$

$$\begin{split}K(T)^{*}&=argmin_{K\in R^{+}}\mathbb{E}\Big{[}(P_{o}+P_{u})\max\{S(T)-K,0\}\\ &+P_{o}(K-S(T))|\mathcal{F}_{t}\Big{]},\end{split}$$

$$\begin{split}\displaystyle K(T)^{*}&=argmin_{K\in R^{+}}\mathbb{E}\Big{[}P_{o}\Big{(}\frac{P_{o}+P_{u}}{P_{o}}\max\{S(T)-K,0\}\\ &+(K-S(T))\Big{)}|\mathcal{F}_{t}\Big{]}.\end{split}$$

$$\mathcal{X}=\Phi (S(T)),$$

$$\Phi (S(T))=\max \{S(T)-K,0\}$$

$$\Pi (t;\mathcal{X})=F(t,S(t)).$$

$$\begin{array}{rl}& \frac{\partial F(t,S(t))}{\partial t}+\mu (t)S(t)\frac{\partial F(t,S(t))}{\partial S}\\ & +\frac{1}{2}S(t{)}^2{\sigma }^2(t))\frac{\partial F(t,S(t))}{\partial S^2}-rF(t,S(t))=0\end{array}$$

$$F(T,S(T))$$

$$\displaystyle F(t,S(t))=e^{-r\Delta t}\mathbb{E}_{S}\big{[}\Phi(S(T))\ \mid\mathcal{F}_{t}\big{]}.$$

$$f(y)=\frac{1}{\sigma (t)\sqrt{2\pi \Delta t}}{e}^{-}\left(\frac{(y-(\mu (t)-\frac{1}{2}{\sigma }^2(t))\Delta t{)}^2}{2\sigma (t{)}^2\Delta t}\right).$$

$$\displaystyle\mathbb{E}_{S}\big{[}\Phi(S(T))\ \mid\mathcal{F}_{t}\big{]}$$

$$\displaystyle=\mathbb{E}_{S}\big{[}\max\{S(T)-K,0\}|\mathcal{F}_{t}\big{]}$$

$$\displaystyle=\mathbb{E}_{S}\big{[}\max\{S(t)e^{y}-K,0\}\big{]}$$

$$={\int }_{ln\frac{K}{S(t)}}^{\infty }S(t)e^{y}f(y)dy-{\int }_{ln\frac{K}{S(t)}}^{\infty }Kf(y)dy$$

$$I_2=$$

$$=$$

$$I_1=$$

$$=$$

$$F(t,S(t))=e^{\mu (t)\Delta t}\mathcal{N}(d_1)S(t)-\mathcal{N}(d_2)K.$$

$$\omega =\frac{P_u}{P_o}$$

$$K^{\ast }(T)=$$

$$=$$

$$\begin{split}\ln\left(\frac{S(t+\Delta t)}{S(t)}\right)&=\left(\mu(t)-\frac{1}{2}\sigma^{2}(t)\right)\Delta t\\ &+\sigma(t)\big{(}W(t+\Delta t)-W(t)\big{)}.\end{split}$$

$$\begin{split}\hat{\sigma}(t)=&\Bigg{(}\frac{1}{n}\sum_{i=2}^{n}\Big{(}\ln\Big{(}\frac{S(t-(n-i)\Delta t)}{S(t-(n-i+1)\Delta t)}\Big{)}\\ -&\frac{1}{n}\sum_{i=1}^{n-1}\Big{[}\ln\Big{(}\frac{S(t-(n-i)\Delta t)}{S(t-(n-i+1)\Delta t)}\Big{)}\Big{]}\Big{)}^{2}\Bigg{)}^{\frac{1}{2}},\end{split}$$

$$\begin{split}\hat{\mu}(t)&=\frac{1}{n}\sum_{i=1}^{n}\Big{[}\ln\Big{(}\frac{S(t-(n-i)\Delta t}{S(t-(n-i+1)\Delta t)}\Big{)}\Big{]}\\ &+\frac{1}{2}\hat{\sigma}(t)^{2},\end{split}$$

$$\begin{split}\textrm{WMAE}&=\frac{1}{N}\sum_{t=1}^{N}\Big{(}\mathbbm{1}_{(S(t)>K^{*}(t))}\omega|S(t)-K^{*}(t)|\\ &+\mathbbm{1}_{(S(t)<K^{*}(t)}|S(t)-K^{*}(t)|\Big{)}\end{split}$$

$$\begin{split}\textrm{WMAPE}&=\frac{1}{N}\sum_{t=1}^{N}\Big{(}\frac{\mathbbm{1}_{(S(t)>K^{*}(t))}\omega|S(t)-K^{*}(t)|}{S(t)}\\ &+\frac{\mathbbm{1}_{(S(t)<K^{*}(t))}|S(t)-K^{*}(t)|}{S(t)}\Big{)},\end{split}$$

Full text

I Introduction

In many applications including manufacturing, energy, and finance, accurate prediction is required to support strategic, tactical and/or operational decisions of organization [1]. When physical information about the underlying mechanism that generates the time series data is limited, data-driven methods can be useful for predicting future observations [2]. In general, data-driven forecasting methods predict future observations based on past observations [3]. Several data-driven methods have been proposed in the literature for modeling time series data, among which Auto-Regressive Integrated Moving Average (ARIMA) and its variants such as the ARIMA-General Auto Regressive Conditional Heteroskedasticity (ARIMA-GARCH) have been widely used in many applications due to their flexibility and statistical properties [4, 5, 6, 7]. ARIMA which assumes a constant standard deviation of stochastic noises, whereas ARIMA-GARCH extends it by allowing the standard deviation to vary over time.

The typical ARIMA-based models estimate its model parameters using historical data and uses the estimated time-invariant parameters throughout the prediction period. Using such time-invariant parameters may not capture possible changes in the underlying data generation mechanism. Some studies modify the original ARIMA model to update the parameters using new observations [8, 9]. The basic idea of these ARIMA-based models is that the future observation can be predicted by using a linear combination of past observations (and estimated noises). Therefore they assume a linear correlation structure between consecutive observations [10]. However, when the underlying dynamics exhibits a highly volatile process, such a simple linear structure may provide poor prediction performance [11].

This study aims to provide accurate predictions for a highly volatile and time-varying stochastic process whose underlying dynamics is complicated and possibly nonlinear. As an example, let us consider a prediction problem faced by a contract manufacturer (CM) located in Michigan in the U.S, which motivates this study. The CM is a manufacturing company that produces various automotive parts, such as front and rear bumper beams, for several large automotive companies worldwide. The CM deals with a large number of orders for bumper beams from several automotive companies and the order sizes are time-varying. The CM should plan their production capacity carefully so that it can deliver products promptly when it gets orders. When an actual order size is greater than expected (i.e., when an order size is underestimated), overtime wages must be paid to workers to meet demands. On the other hand, when an order size is smaller than predicted (i.e., when an order size is overestimated), workers and equipment become idle.

As such, CM wants to predict future order sizes accurately, so that it can reduce its operating costs resulting from the discrepancy between its predicted value and actual sizes. Currently, CM uses its own proprietary prediction model, but its prediction performance is not satisfactory. The details of CM’s proprietary model are confidential, so we cannot find reasons for its unsatisfactory performance. When we apply the ARIMA and ARIMA-GARCH models to CM’s datasets, we also do not obtain significantly better prediction results (detailed results will be provided in Section  III). We believe such poor performance of ARIMA-based approaches is because they cannot fully characterize the underlying volatile dynamics. In addition to historical data, the future order size may depend on other factors which possibly make the order process behave nonlinearly. A new prediction approach that can adapt to such time-varying, and possibly nonlinear, dynamics is needed for providing better forecasts.

To this end we develop a new method for predicting future values in highly volatile processes, based on real option pricing theories typically used in financial engineering. One of the popularly used stochastic process models for pricing real options is the Geometric Brownian Motion (GBM) model. Brownian motion is a continuous-time stochastic process, describing random movements in time series variables. The GBM, which is a stochastic differential equation, incorporates the idea of Brownian motion and consists of two terms: a deterministic term to characterize the main trend over time and a stochastic term to account for random variations. In GBM the random variations are represented by Brownian Motion [12]. GBM is useful to model a positive quantity whose changes over equal and non-overlapping time intervals are identically distributed and independent.

The GBM has been applied to represent various real processes in finance, physics, etc. [13]. In particular, it becomes a fundamental block for many asset pricing models [12], and recently it has been applied to facilitate the use of a rich area of options theory to solve various pricing problems (see, for example, [14, 15, 16, 17, 18, 19]). However, most of the current GBM studies have been limited to solving pricing problems and have not used real options theory for making forecasts.

In this study, by utilizing the full power of real options theory, we present a new approach for predicting future observations when the system’s underlying dynamics follows the GBM process. Specifically, we allow the GBM parameters to adaptively change over time in order to characterize time-varying dynamics. We formulate the prediction problem as an optimization problem and provide a solution using real option theories. To the best of our knowledge, our study is the first attempt to incorporate options theory in the prediction problem.

Our approach provides extra flexibility by allowing overestimation (or over-prediction) to be handled differently from underestimation (or over-prediction). The overestimation and underestimation costs are determined in real life applications, depending on a decision-maker’s (or organization’s) preference. For example, in the aforementioned CM case, overestimation and underestimation of order sizes could cause different costs. The CM may want to put a larger penalty on the demand underestimation than on the overestimation, so that it can avoid extra overtime wages. We incorporate unequal overestimation and underestimation costs into the optimization problem and find the optimal forecast that minimizes the expected prediction cost.

To evaluate the prediction performance, we use three datasets collected from different applications, including the demand for bumper beams in CM (manufacturing), stock prices (finance), and wind speed (environment). We compare the performance of our model with ARIMA and ARIMA-GARCH models (and the proprietary prediction model in the CM case study) with different combinations of overestimation and underestimation costs. In most cases, our model outperforms those alternative models. In particular, we find that when the process is highly time-varying such as stock prices and wind speed, the proposed approach provides much stronger prediction capability than ARMA and ARIMA-GARCH.

The remainder of the paper is organized as follows. The mathematical formulation and solution procedure are discussed in Section II. Section III provides numerical results in three different applications. Section IV concludes the paper.

II Methodology

II-A Problem Formulation

Consider a real-valued variable $S(t)$ which represents a system state at time $t$. For example, the state variable can be a stock market index price, a manufacturer’s order size, or wind speed. This state variable is assumed to follow an inhomogeneous GBM with time-varying parameters.

Let us consider a filtered probability space $(\Omega,\mathcal{F},P,\mathcal{F}_{t})$, where the filtration $\mathcal{F}_{t}$ is generated by the Brownian motion $W$, i.e. $\mathcal{F}_{t}=\mathcal{F}_{t}^{{W}}$ so that $\mathcal{F}_{t}$ contains all information generated by ${W(t)}$, up to and including time $t$. With GBM, the stochastic process $S(t)$ is modeled by the following dynamics.

[TABLE]

where $\sigma(t)$ denotes the volatility of $S(t)$ and $\mu(t)$ represents a drift process. The stochastic process $W(t)$ represents the Brownian motion where the increment $W(t+\Delta t)-W(t)$ during the time interval $\Delta t$ is normally distributed with mean 0 and variance $\Delta t$, denoted by $\mathcal{N}(0,\Delta t)$, and $W(t)$ is assumed to be stationary.

Our objective is to predict $S(T)$ in the future time at $T(>t)$ when the current time is $t$. Solving (1) by using Itô’s lemma [20], we obtain

[TABLE]

and

[TABLE]

Let $K$ be the predicted value of $S(T)$ at time $T$. When the overestimation and underestimation is penalized equally, the quantity that represent the variable’s central tendency, such as mean and median, is commonly used for prediction. But we consider a more general case where overestimation needs to penalized differently from underestimation, as discussed in Section I. When the observed value is $S(T)$, the overestimated quantity becomes $\max\{K-S(T),0\}$, while the underestimated quantity is $\max\{S(T)-K,0\}$.

Let $p_{o}$ and $p_{u}$ denote the penalties for over/underestimation, respectively. We formulate the optimization problem for estimating $S(K)$ that can minimize the expected prediction cost,

[TABLE]

Note that

[TABLE]

If we substitute (5) into (4), the optimal predicted value, denoted by $K^{*}$, can be obtained by solving the following objective function.

[TABLE]

or equivalently,

[TABLE]

In the next section we will present a solution procedure to obtain $K^{*}(T)$, based on the option theory.

II-B Real Option Based Solution Procedure

The optimization problem in (7) can be reformulated by employing the financial pricing theories. Suppose that we want to predict a state at the future time $T$. In pricing theories, $T$ can be viewed as the date to maturity, or the expiration date.

A real option, also called contingent claim, with the date to maturity $T$, can be constructed on the state variable $S(t)$. A real option is a stochastic variable $\mathcal{X}\in\mathcal{F}_{T}^{W}$ that can be expressed as

[TABLE]

where $\Phi(\cdot)$ is a contract function.

The contract function $\Phi(\cdot)$ is typically set to the payoff of the real option at time $T$. When the predicted value is $K$, $K$ can be viewed as the strike value in the option theory, while $\max\{S(T)-K,0\}$ is the payoff. Therefore, we get

[TABLE]

It is required that $\mathcal{X}\in\mathcal{F}_{T}^{W}$ ensures that the value of the payoff of the real option $\mathcal{X}$ is determined at time $T$.

Let the price process $\Pi(t;\mathcal{X})$ for the real option at time $t$ be given by a function $F(t,S(t))\in[t,T]\times R_{+}$, i.e.,

[TABLE]

Here $F(\cdot)$ is a function which is assumed to be once continuously differentiable in $t$, and twice in $S(t)$.

For a short-term prediction, the time interval $\Delta t$ between the current time $t$ and the future time $T$ is small, so we can assume that $\mu(t)$ and $\sigma(t)$ are constants during $[t,T]$. Then $F(t,S(t))$ can be obtained by solving the Black-Scholes Partial Differential Equation (PDE) [20],

[TABLE]

with

[TABLE]

where $r$ represents a discounting factor.

The Black-Scholes PDE in (11)-(12) is usually solved numerically. But alternatively, we solve it using the Feyman-Ka$\breve{c}$ stochastic representation formula [20], to obtain

[TABLE]

Next, we derive $F(t,S(t))$ in a closed form, given $K$. Letting y=ln\Big{[}{S(T)}/{S(t)}\Big{]} and using the fact that $S(T)=S(t)exp((\mu(t)-\frac{1}{2}\sigma^{2}(t))\Delta t+\sigma(t)\Delta W(t))$, it follows that y\sim\mathcal{N}\big{(}\big{(}\mu(t)-\frac{1}{2}\sigma^{2}(t)\big{)}\Delta t,\sigma^{2}(t)\Delta t\big{)}. Thus, the probability density function $f(y)$ of $y$ is given by

[TABLE]

Consequently, we obtain

[TABLE]

To solve (17), let $I_{1}$ and $I_{2}$, respectively, denote the first and second terms in (17). We also let $z={y-(\mu(t)-0.5\sigma^{2}(t))\Delta t}/{\sigma(t)\sqrt{\Delta t}}$. First, $I_{2}$ becomes

[TABLE]

where $d_{2}={ln(\frac{S(t)}{K})+(\mu(t)-\frac{1}{2}\sigma^{2})\Delta t}/{\sigma(t)\sqrt{\Delta t}}$ and $\mathcal{N}(\cdot)$ denotes the cumulative distribution function (CFD) for the standard normal distribution. Next, we obtain $I_{1}$ as

[TABLE]

where we use $v=z-\sigma(t)\sqrt{\Delta t}$ in (24) and $d_{1}=d_{2}+\sigma(t)\sqrt{\Delta t}$ in (25).

For small $\Delta t$, we can set $r=0$. Then, $F(t,S(t))$ in (13) becomes:

[TABLE]

Note that given $\mu(t)$ and $S(t)$ at the current time $t$ and $K$, we can obtain $F(t,S(t))$.

With the obtained expected payoff \mathbb{E}_{S}\big{[}\Phi(S(T))\ \mid\mathcal{F}_{t}\big{]} where $\Phi(S(T))=\max\{S(T)-K,0\}$, we can find the optimal $K^{*}(T)$ in (7). Let $\omega$ denote the ratio of overestimation cost to underestimation cost, i.e.,

[TABLE]

Given the price of the real option, defined in (13), we can reformulate the optimization problem in (7) as

{strip}

[TABLE]

with $d_{2}={ln(\frac{S(t)}{K})+(\mu(t)-\frac{1}{2}\sigma^{2})\Delta t}/{\sigma(t)\sqrt{\Delta t}}$ and $d_{1}=d_{2}+\sigma(t)\sqrt{\Delta t}$. We use (13) with $r=0$ in the first term in the second equality and the last term in the second equality is obtained using (3). By plugging $F(t,S(t))$ in (26), we get the last equality.

The predictor $K^{*}$ prefers overestimation when $\omega>1$ or underestimation when $\omega<1$. When overestimation and underestimation are equally penalized, the optimal $K^{*}$ can be obtained with $w=1$ in (28). The optimization function in (28) is a convex optimization problem that can be solved efficiently by existing numerical optimization softwares. In our implementation, we use Scipy’s (Scientific Python) optimization library in Python.

II-C Parameters Estimation

For a volatile stochastic process, the parameters $\mu(t)$ and $\sigma(t)$ can be time-varying. We estimate the nonstationary parameters using recent observations. Consider $n$ recent observations at the current time $t$, i.e., $S(t-(n-1)\Delta t),S(t-(n-2)\Delta t),\cdots,S(t)$. Because $S(t)$ follows geometric Brownian motion and $\mu(t)$ and $\sigma(t)$ are assumed to be constant during the short interval $\Delta t$, the discretization scheme of (2) is given by

[TABLE]

Noting that under GBM \ln\Big{(}{S(t+\Delta t)}/{S(t)}\Big{)} is normally distributed with mean \big{[}\mu(t)-\frac{1}{2}\sigma^{2}(t)\big{]}\Delta t and variance $\sigma^{2}(t)$, we estimate $\mu(t)$ and $\sigma(t)$ using maximum likelihood method as

[TABLE]

[TABLE]

respectively.

The estimated parameters $\hat{\mu}(t)$ and $\hat{\sigma}(t)$ are plugged into (26) and we obtain the optimal predicted value $K^{*}$ for $S(t+\Delta t)$ by solving (30).

II-D Implementation Details

We refer our proposed model to as the option prediction model. Figure 1 summarizes the overall procedure of the proposed approach. We also summarize the procedure of the proposed approach in Algorithm 1 below. We set the time step $\Delta t=1$ to make the one-ahead step prediction. The data is divided into three sets: training, validation, and testing. The training set starts at $t=1$ and ends at $t=N_{1}$, consisting of about 50% of the entire data set, is used to determine the model parameters as shown in Figure 1. The validation set, consisting of about 20% of the data set, is used for determining the window size $n$. Lastly the testing set consists of the last $30\%$ of the data set and it starts at $t=N_{2}$.

In Algorithm 1 we determine the window size $n$ for obtaining the parameters $\hat{\mu}(t)$ and $\hat{\sigma}(t)$, we use the validation technique [21]. We fit the model with a different window size $n$ and evaluate the prediction performance using data in the validation set and choose the best window size that generates the lowest prediction error in the validation set. The performance of our approach is evaluated using data in the testing set (See Figure 1). We report the prediction performance in the testing set in Section III.

In evaluating the prediction performance, we consider that the overestimated and underestimated prediction results need to be evaluated differently for $\omega\neq 1$. As such we employ the following two performance measures, namely, Weighted Mean Absolute Error (WMAE) and Weighted Mean Absolute Percentage Error (WMAPE), defined by

[TABLE]

and

[TABLE]

respectively, where $N$ denotes the number of data points in the testing set and $K^{*}(t)$ is the predicted value at time $t$.

III Case Studies

This section implements the proposed prediction model using multiple datasets obtained from real-life applications. Specifically we examine the performance of the predictive model in predicting the size of a manufacturer’s order, a stock market index price, and wind speed.

III-A Alternative methods

We compare our model with two standard time series models, namely, the ARIMA and the ARIMA-GARCH. We use the Akaika Information Criteria (AIC) to select the model order in both models. For fair comparison, we update the model parameters in a rolling horizon manner, similar to the procedure discussed in Section II-D. That is, we determine the window size $n$ using the validation technique and update the model parameters using the most recent $n$ observations whenever a new observation is obtained.

With underestimation penalties, Pourhab et al. [22] suggest using quantile of the predictive state density. With $\omega(=p_{u}/p_{o})$ denoting the ratio of underestimation cost to overestimation cost, we use the $({\omega}/{1+\omega})$-quantile, given by

[TABLE]

where $\hat{\mu}_{a}(t)$ denotes the estimated predicted mean, $\hat{\sigma}_{a}$ is the estimated standard deviation in ARIMA (or ARIMA-GARCH) model, and $\Phi^{-1}(\cdot)$ denotes the inverse of the standard normal CDF. Note that large (small) $w$ puts more penalty on $p_{u}$ ($p_{o}$) and the quantile prediction provides a larger (smaller) predicted value, so underestimation (overestimation) can be avoided.

III-B Manufacturing Data

We first study the prediction problem faced by our industry partner, CM. The historical data obtained from CM includes orders of 10 different types of bumper beams. We use monthly data on those 10 types of bumper beams ordered over a period of 29 consecutive months (the order size varies from 0 to over 36,000 items). When applying the proposed model to this problem, the choice of weight $\omega$ affects the final prediction. By changing the weight, we are able to show a preference for over-capacity (overestimation) or under-capacity (underestimation). We consider different cases for choosing the weight parameter $\omega$.

Let us first look at the case when $\omega$ is set to be less than one (i.e, $p_{u}\leq p_{o}$). According to CM, workers and equipment can be shifted from one type of bumper beam to another, but doing so incurs 15$\%$ loss of production efficiency. In other words, if one type of bumper beam is overestimated, causing over-capacity, available resources can be assigned to other bumper beam production, but with a reduced efficiency. In this case, underestimation is favored and we set $w=1/1.15$.

Next the weight parameter can be set to be greater than 1 (i.e, $p_{u}\geq p_{o}$) when the prediction is preferred to be more than the actual order size. According to the labor law in Michigan in the U.S., overtime rate is higher than the regular salary. In this case we set $w=1.15$ to emphasize the preference of overestimation to underestimation. Finally we also consider $w=1$, which reflects equal penalties.

The errors in terms of WMAE and WMAPE for all ten types of bumper beams are presented in Tables I-III with three different weights. Overall our option prediction model performs better than the CM’s own prediction, ARIMA and ARIMA-GARCH in both criteria. With $w=1/1.15$ the proposed approach provides lower WMAEs (WMAPEs) for 9 (5) types of bumper beams out of 10 types. Similarly, with other $w$ values, our approach outperforms the alternative models in most cases.

Although ARIMA and ARIMA-GARCH provide the lowest errors for some products, their prediction performance is not consistent. For example, for $1^{st}$, $4^{th}$ and $7^{th}$ product, WMAEs from ARMA are much higher than the proposed approach, whereas ARIMA-GARCH results in pretty poor performance for predicting order sizes for $8^{th}-10^{th}$ products. On the contrary, our approach provides more stable results. Even when WMAEs and WMAPEs from our approach are higher than other approaches, they are close to the lowest errors. Therefore, we can conclude that our approach is more accurate and reliable. The CM’s proprietary model does not account for unequal weights on overestimation and underestimation. If the company wants to minimize the excess inventory due to overestimation, a small (less than 1) weight parameter should be assigned. If the company goal is to meet customer satisfaction, overestimation should be preferred with a large (larger than 1) weight parameter. In this sense our approach can reflect the company’s management preference more flexibly.

III-C Stock Market Index Data

To evaluate the performance of our approach in a highly volatile process, we consider stock market index price time series data. We analyze the daily closing price of the Dow Jones index in three time periods between 2010 and 2015.

Risk averse and risk seeking investors have different preferences in terms of overestimation and underestimation. That being said, in a bull market, stock prices are expected to increase. In such a case, risk seeking investors with aggressive investment strategies would prefer biasing their prediction to overestimation. On the contrary, risk averse investors tend to be less optimistic, making them conservative, preferring underestimation. To reflect different investment preferences, we consider three values of the weight parameter $\omega$, $1/1.15$, $1$, or $1.15$, to represent the underestimation preference, neutral/no preference, and overestimation preference, respectively.

Table IV summarizes the results with three testing periods. Each testing period includes 100 days. Clearly, our option prediction performs better than ARIMA-GARCH and ARIMA in all cases, alerting for the possibility of a profitable trading strategy. The ARMA and ARMA-GARCH models generate 2.5 to 10 times higher WMAEs and 2 to 11 times higher WMAPEs.

III-D Wind Speed Data

Finally, we additionally consider another highly volatile process, wind speed. Because of environmental considerations, wind power, as a renewable source of energy, has been increasingly adopted worldwide [23]. Intermittent output of the farm is considered a challenging issue in terms of integrating the wind power into electric power grids. For reliable supply of power, steady and uninterrupted energy generation is desirable, which is not the case with wind energy. Wind speed is highly variable, depending on weather conditions and geographical factors such as the terrain. Such variability imposes challenges in power grid operations. To overcome the challenges, accurate forecasting of wind speed is required [24].

We use wind speed data collected from a meteorological tower near a wind farm located in Europe. The whole dataset consists of about 3000 samples, which covers a period of about a month. Due to the data confidentiality required by our industry partner, we omit more detailed description of the dataset studied in this case study. In wind farm operations some operators want to put a higher penalty on overestimation to avoid unsatisfied demand (or unsatisfied commitment), whereas underestimating wind speeds may be preferred when the salvage cost of excessively generated power is high [22, 25]. To reflect different costs, we use three different values for $\omega$, $1/1.15$, $1$ and $1.15$.

Table V summarizes the prediction results in the testing set from the three models. The proposed option prediction significantly outperform the other methods. The WMAEs and WMAPEs from ARMA and ARIMA-GARCH are higher by one order of magnitude than our approach. It demonstrates the superior prediction performance of our approach in a highly volatile process.

IV Conclusion

In this study, we present a new prediction methodology for the time series data, based on option theories in finance when the underlying dynamics is assumed to follow the GBM process. To characterize time-varying patterns, we allow the GBM model parameters to vary over time and update the parameter values using recent observations. We formulate the prediction problem with unequal overestimation and underestimation penalties as the stochastic optimization problem and provide its solution procedure. We demonstrate the prediction capability of the proposed approach in various applications. Our approach appears to work well in the manufacturing application, when the order size varies over time. For more highly volatile processes such as stock prices and wind speeds, the proposed model exhibits much stronger prediction capability, compared to alternative ARIMA-based models.

In the future, we plan to investigate other parameter updating schemes. In this study, we update parameters in a rolling horizon manner using the maximum likelihood estimations. Another possibility is to use the Kalman filtering or its variants. Long-term predictions are beyond the scope of this study, but we plan to extend the approach presented in this study for obtaining accurate long-term predictions. We will also incorporate prediction results into managerial decision-making in several applications such as power grid operation with renewable energy [26].