Regularized two-stage stochastic variational inequalities for Cournot-Nash equilibrium under uncertainty

TL;DR

This paper formulates a two-stage stochastic variational inequality model for Cournot-Nash equilibrium under uncertainty, providing solution existence conditions, a regularized approximation method, and demonstrating its application to oil market data.

Contribution

It introduces a regularized sample average approximation method for two-stage stochastic variational inequalities and proves its convergence, applied to oligopolistic market modeling.

Findings

Proven convergence of the proposed method as regularization and sample size increase.

Effective modeling of oil market shares using the two-stage SVI approach.

Numerical validation with historical crude oil data.

Abstract

A convex two-stage non-cooperative multi-agent game under uncertainty is formulated as a two-stage stochastic variational inequality (SVI). Under standard assumptions, we provide sufficient conditions for the existence of solutions of the two-stage SVI and propose a regularized sample average approximation method for solving it. We prove the convergence of the method as the regularization parameter tends to zero and the sample size tends to infinity. Moreover, our approach is applied to a two-stage stochastic production and supply planning problem with homogeneous commodity in an oligopolistic market. Numerical results based on historical data in crude oil market are presented to demonstrate the effectiveness of the two-stage SVI in describing the market share of oil producing agents.

Tables

Table 1. Table 1: Numerical results for different ϵ italic-ϵ \epsilon and sample size ν 𝜈 \nu , J = 10 𝐽 10 J=10 with individual sample.

$\nu$	$J\left(1+2\nu \right)$	Iter	CPU time/s	$\mathrm{𝐑𝐞𝐬}$	Iter	CPU time/s	$\mathrm{𝐑𝐞𝐬}$
		$ϵ={10}^{-3}$			$ϵ={10}^{-6}$
10	210	146.30	0.26	4.42e-01	176.20	0.32	3.88e-04
50	1010	194.70	1.81	9.35e-01	197.40	1.83	9.32e-04
500	10010	208.70	26.72	3.00e+00	212.20	27.21	2.99e-03
2000	40010	222.60	154.97	5.93e+00	220.50	153.54	6.00e-03
5000	100010	224.70	623.53	9.49e+00	226.40	627.53	9.48e-03
		$ϵ={10}^{-9}$			$ϵ={10}^{-12}$
10	210	152.70	0.27	1.08e-06	169.40	0.30	9.49e-07
50	1010	197.20	1.83	1.41e-06	194.40	1.80	9.75e-07
500	10010	212.70	27.21	3.21e-06	209.70	26.85	9.59e-07
2000	40010	220.30	153.34	6.16e-06	220.70	153.73	9.51e-07
5000	100010	226.70	628.89	9.60e-06	226.20	627.58	9.60e-07

Table 2. Table 2: Numerical results for different ϵ italic-ϵ \epsilon and sample size ν 𝜈 \nu , J = 10 𝐽 10 J=10 with block implementation.

$\nu$	$J\left(1+2\nu \right)$	Iter	CPU time/s	$\mathrm{𝐑𝐞𝐬}$	Iter	CPU time/s	$\mathrm{𝐑𝐞𝐬}$
		$ϵ={10}^{-3}$			$ϵ={10}^{-6}$
10	210	162.60	0.21	4.04e-01	155.50	0.20	4.14e-04
50	1010	180.60	0.94	9.42e-01	184.00	0.93	9.31e-04
500	10010	203.30	10.67	2.97e+00	204.10	10.61	2.98e-03
2000	40010	213.80	42.38	5.89e+00	213.70	42.36	5.91e-03
5000	100010	218.20	115.15	9.36e+00	219.60	115.63	9.36e-03
		$ϵ={10}^{-9}$			$ϵ={10}^{-12}$
10	210	143.20	0.19	1.08e-06	149.50	0.19	9.59e-07
50	1010	197.40	1.02	1.38e-06	191.10	0.97	9.58e-07
500	10010	205.10	10.58	3.19e-06	202.40	10.67	9.52e-07
2000	40010	214.10	42.46	6.07e-06	213.40	42.62	9.73e-07
5000	100010	219.00	115.35	9.50e-06	219.70	115.58	9.74e-07

Equations

$$(x_j\in \mathbb{R},y_j:\Xi \to \mathbb{R},\mathcal{F}\text{-measurable}),$$

$$(x\in {\mathbb{R}}^J,y:\Xi \to {\mathbb{R}}^J,\mathcal{F}\text{-measurable}).$$

$${\text{maximize}}_{(x_j,y_j(\cdot ))}$$

$$x_j\in X_j,$$

$$\displaystyle y_{j}(\xi)\in_{a.s.}Y_{j},~{}~{}g_{j}\big{(}\xi,x_{j},y_{j}(\xi)\big{)}\leq_{a.s.}0$$

$$x_{-j}^{\ast }$$

$$y_{-j}^{\ast }(\xi )$$

$$\text{maximize}~{}W^{1}_{j}(x_{j},x^{*}_{-j})+\mathbb{E}\big{[}\Phi_{j}(x_{j},x^{*}_{-j},\xi)\big{]},$$

$${\Phi }_j(x_j,x_{-j}^{\ast },\xi )=\underset{y_j(\xi )\ge 0}{\sup }\{W_j^2(\xi ,y_j(\xi ),y_{-j}^{\ast }(\xi ))|x_j\ge y_j(\xi ),\text{ for a.e. }\xi \in \Xi \}.$$

$$L_{j}(x_{j},x^{*}_{-j},y_{j},y^{*}_{-j},\lambda_{j})=L_{j}^{1}(x_{j},x^{*}_{-j})+\mathbb{E}\big{[}L_{j}^{2}\big{(}\xi,x_{j},y_{j}(\xi),y^{*}_{-j}(\xi),\lambda_{j}(\xi)\big{)}\big{]},$$

$$L_j^1(x_j,x_{-j}^{\ast })=W_j^1(x_j,x_{-j}^{\ast }),$$

$$\displaystyle L_{j}^{2}\big{(}\xi,x_{j},y_{j}(\xi),y^{*}_{-j}(\xi),\lambda_{j}(\xi)\big{)}=W^{2}_{j}(\xi,y_{j}(\xi),y^{*}_{-j}(\xi))+\lambda_{j}(\xi)(x_{j}-y_{j}(\xi)).$$

$$W_j^1(x_j,x_{-j}^{\ast })=-\frac{1}{2}{c}_j{x}_j^2-a_j{x}_j,$$

$$W^{2}_{j}(\xi,y_{j}(\xi),y_{-j}(\xi))=\big{(}p_{j}(\xi)-\gamma(\xi)T(y(\xi))\big{)}y_{j}(\xi),$$

$${\text{maximize}}_{x_j}$$

$$0\le x_j,$$

$${\Phi }_j(\xi ,x)=$$

$$0\le y_j(\xi )\le x_j,\text{for a.e. ξ∈Ξ}.$$

$$\displaystyle-c_{j}\bar{x}_{j}-a_{j}+\mathbb{E}\big{[}\bar{\lambda}_{j}(\xi)\big{]}\in\mathcal{N}_{[0,\infty)}(\bar{x}_{j}),$$

$$p_j(\xi )-\gamma (\xi )\sum _{i\ne j}^J{\bar{y}}_i(\xi )-2\gamma (\xi ){\bar{y}}_j(\xi )-{\bar{\lambda }}_j(\xi ){\in }_{a.s.}{\mathcal{N}}_{[0,\infty )}({\bar{y}}_j(\xi )),$$

$${\bar{x}}_j\ge 0,{\bar{y}}_j(\xi ){\ge }_{a.s.}0,{\bar{x}}_j-{\bar{y}}_j(\xi ){\ge }_{a.s.}0,$$

$${\bar{\lambda }}_j(\xi ){\ge }_{a.s.}0,$$

$${\bar{\lambda }}_j(\xi ){\perp }_{a.s.}({\bar{x}}_j-{\bar{y}}_j(\xi )).$$

$$0\le x_j$$

$$0{\le }_{a.s.}{y}_j(\xi )$$

$$0{\le }_{a.s.}{\lambda }_j(\xi )$$

$$0$$

$$C=\text{diag}(c_1,c_2,...,c_J),\Pi (\xi )=\gamma (\xi )(e{e}^T+I).$$

$$\big{(}x^{*},y^{*},\lambda^{*}\big{)}=\big{(}(x_{1}^{*},y_{1}^{*},\lambda_{1}^{*}),\ldots,(x_{J}^{*},y_{J}^{*},\lambda_{J}^{*})\big{)}\in\mathbb{R}^{J}\times\mathfrak{L}^{\infty}_{J}\times\mathfrak{L}^{1}_{J}$$

$$M(\xi )=\left(\begin{array}{cc}\Pi (\xi )& I\\ -I& 0\end{array}\right)$$

$$\left(\begin{array}{c}\hat{y}(\xi )\\ \hat{\lambda }(\xi )\end{array}\right)\ge 0,\left(\begin{array}{cc}\Pi (\xi )& I\\ -I& 0\end{array}\right)\left(\begin{array}{c}\hat{y}(\xi )\\ \hat{\lambda }(x)\end{array}\right)+\left(\begin{array}{c}-p(\xi )\\ x\end{array}\right)\ge 0,\text{ for a.e. ξ∈Ξ}.$$

$$0\le \left(\begin{array}{c}{y}_1(\xi )\\ y_2(\xi )\\ {\lambda }_1(\xi )\\ {\lambda }_2(\xi )\end{array}\right)\perp \left(\begin{array}{cccc}2\gamma (\xi )& \gamma (\xi )& 1& 0\\ \gamma (\xi )& 2\gamma (\xi )& 0& 1\\ -1& 0& 0& 0\\ 0& -1& 0& 0\end{array}\right)\left(\begin{array}{c}{y}_1(\xi )\\ y_2(\xi )\\ {\lambda }_1(\xi )\\ {\lambda }_2(\xi )\end{array}\right)+\left(\begin{array}{c}-p_1(\xi )\\ -p_2(\xi )\\ x_1\\ x_2\end{array}\right)\ge 0,\text{for a.e. ξ∈Ξ}.$$

$$\left\{(0,0,{\tilde{\lambda }}_1(\xi ),{\tilde{\lambda }}_2(\xi )):{\tilde{\lambda }}_1(\xi )=\{\begin{array}{ll}0,& x_1>0\\ {\lambda }_1(\xi ),& x_1=0\end{array},{\tilde{\lambda }}_2(\xi )=\{\begin{array}{ll}0,& x_2>0\\ {\lambda }_2(\xi ),& x_2=0\end{array},\text{for a.e. ξ∈Ξ}\right\},$$

$$M^{ϵ}(\xi )=\left(\begin{array}{cc}\Pi (\xi )& I\\ -I& ϵI\end{array}\right)\mathrm{and}q(x,\xi )=\left(\begin{array}{c}-p(\xi )\\ x\end{array}\right).$$

$$0$$

Full text

0C15, 90C33

1 Introduction

We consider a two-stage non-cooperative multi-agent game under uncertainty. The model is used to describe a homogenous commodity production and supply planning problem in an oligopolistic market. In particular, we focus on a $J$-player stochastic Nash equilibrium problem (SNEP) of Cournot competition, whose solution concept is characterized by stochastic Cournot-Nash (CN) equilibria. Conventional non-cooperative game theory has a long history being an effective model to describe market behaviour, yet mostly under deterministic settings. In order to explore characteristics of real markets, one cannot neglect the presence of uncertainty.

Researchers have studied various real markets through SNEPs in recent years. Jofré, Rockafellar and Wets [24] investigated various economic equilibria using SVIs. A scenario-based multi-stage oligopolistic market equilibrium problem under uncertainty was discussed in [16]. A two-settlement oligopolistic equilibrium with uncertainty in the future market was presented in [48]. For practical applications, electricity markets with hydro-electric distribution have been studied by Philpott, Ferris and Wets [30] in which the levels of water reserves were modelled under uncertainty. The contemporary treatment of classical equilibrium problems are investigated through finite-dimensional variational inequalities (VIs) and complementarity problems (CPs) with a wide range of applications under the assumption of deterministic and single-stage decision, see [20, 14] and references therein. Yet, many practical applications are formulated under uncertainty [21, 10, 29, 47, 9, 28], and the corresponding multi-stage SVIs and CPs were studied extensively over the last two decades. In [5], Chen and Fukushima considered the stochastic linear CP (LCP) by expected residual minimization (ERM) procedure. The quasi-Monte Carlo method was adopted to generate scenarios of observation and thus to obtain the discrete approximation problem. Chen, Wets and Zhang [10] investigated SVI problems by the ERM procedure, and the sample average approximation (SAA) method was employed to approximate the expected smoothing residual function. More recently, as an extension from single-stage to multi-stage decision processes, Rockafellar and Wets [40] put forward multi-stage SVI problems, which laid a solid foundation for further research. In [6], Chen, Pang and Wets introduced the ERM procedure for two-stage SVI problems and the Douglas-Rachford splitting method was used to present numerical results. Chen, Sun and Xu [9] considered a two-stage stochastic LCP. Structural properties of the problem were studied under the assumption of strong monotonicity, and a discrete scheme was conducted by partition of the support set and the corresponding convergence assertion was established. More generally, Chen, Shapiro and Sun [8] investigated the SAA of two-stage stochastic nonlinear generalized equations, which included two-stage nonlinear SVI problems as special cases. Exponential rate of convergence was derived by using the technique of perturbed partial linearization. From the perspective of the numerical calculation, the equilibrium problems are usually rewritten as a minimization problem, mostly nonsmooth. For this class of problems, the smoothing techniques (see [3]) can be employed so that differentiable methods, e.g., Newton’s method, become applicable in solving the smoothing problem, see for instance [12, 13, 4, 11]. For numerical implementation, Rockafellar and Sun extended the well-known progressive hedging method (PHM) for multi-stage stochastic programming problems to multi-stage SVIs in [33].

The decision vectors of production and supply plan problems, so-called strategies, are distinguished into two categories: (i) those of “here-and-now” type, which do not depend on outcomes of random events in the future, and (ii) those treated as responses. The goal of this paper is to establish a model that describes the market mechanism for Cournot competition under uncertainties. The solution concepts of the model assemble the real strategies adopted by participants of the market. We will provide methodology to solve the proposed model and demonstrating its numerical implementation.

We summarize the main contributions of this paper as follows.

• •

A two-stage stochastic Nash equilibrium problem is proposed to model production and supply competition of a homogenous product under uncertainty in an oligopolistic market.

• •

The model is recast as two-stage SVIs whose solutions characterize a CN equilibrium.

• •

A regularized sample average approximation method is proposed to solve the two-stage SVIs with convergence properties under mild assumptions.

• •

The model is tested numerically for its effectiveness. Moreover, it is used to describe the market share observation in the world market of crude oil. We show that our model is not only able to reproduce historical in-sample market share but also capable of making out-of-sample predictions based on real data sets.

The remaining of the paper is organized as follows. In section 2, the two-stage SNEPs are developed and recast into two-stage SVIs. Section 3 contains structural analysis of two-stage SVIs and the corresponding regularized problems. The convergence theorems of the proposed regularized SAA method are presented in section 4. Numerical results based on randomly generated data and simulation on real market data of crude oil are presented and analyzed in section 6. We conclude our study in the final section.

2 Two-stage stochastic Cournot-Nash game

All stochastic models involve inherently “ordered” components over decision horizons. In particular, the decisions in a strategy may respond to the information that is available only at the present stage, which is known as nonanticipativity. In this paper, we consider a two-stage stochastic CN equilibrium problem, which extends the classical deterministic CN equilibrium problem in [27]. Let $\xi:\Omega\to\mathbb{R}^{d}$ be a random variable with support set $\Xi\subseteq\mathbb{R}^{d}$, and let $(\Xi,\mathcal{F},P)$ denote the induced probability space. A strategy pair of agent $j\in\mathcal{J}:=\{1,...,J\}$ is denoted as

[TABLE]

where $x_{j}$ is a first-stage decision vector and $y_{j}\in\mathcal{Y}$ represents a second-stage response function with $\mathcal{Y}$ being the space of $\mathcal{F}$-measurable functions defined on $\Xi$. Let $\mathfrak{L}_{n}$ be the Lebesgue space of $\mathbb{R}^{n}$-valued functions with $\mathfrak{L}^{\infty}_{n}$ denotes the class of measurable essentially bounded functions. Following a similar treatment as in [6], we further require the second-stage response function of random variable to be essentially bounded, i.e., $y_{j}\in\mathfrak{L}^{\infty}_{1}$. Collectively, the vector of strategy pairs of all agents can be written as

[TABLE]

A strategy pair \big{(}x_{j}^{*},y_{j}^{*}\big{)}\in\mathbb{R}\times\mathfrak{L}^{\infty}_{1} is said to be an equilibrium of our stochastic model if it solves the following problem for all agents $j\in\mathcal{J}$.

[TABLE]

where

[TABLE]

and $y_{j}(\xi)$ denotes the value of response $y_{j}$ to realization $\xi$111For ease the exposition, the same notation $\xi$ is used for both a random variable and its realization in $\mathbb{R}^{d}$ without causing confusion on the context., with

• •

$W^{1}_{j}:\mathbb{R}\times\mathbb{R}^{J-1}\to\mathbb{R}$ a first-stage wealth function of agent $j$, concave and continuously differentiable with respect to $x_{j}$.

• •

$W^{2}_{j}:\Xi\times\mathbb{R}\times\mathbb{R}^{J-1}\to\mathbb{R}$ a second-stage wealth function of agent $j$, concave, well-defined and finite.

• •

$X_{j},Y_{j}$ nonempty, closed and convex subsets of $\mathbb{R}$ and the second-stage constraints hold almost surely ($a.s.$).

• •

$g_{j}:\Xi\times\mathbb{R}\times\mathbb{R}\to\mathbb{R}$ a continuously differentiable function with respect to \big{(}x_{j},y_{j}(\xi)\big{)} for almost every (a.e.) $\xi\in\Xi$ and $\mathcal{F}$-measurable.

In this paper, the model (3) is formulated under the assumption that the uncertainty can be described by a random variable $\xi$ with known distribution. From the perspective of the entire system, the market222A “social planner” is commonly termed in literature of economics, and can be interpreted as one individual who oversees the system and assigns strategy to each agents. “chooses” $x\in\mathbb{R}^{J}$ before a realization $\xi\in\Xi$ is revealed and later “selects” $y(\xi)\in\mathbb{R}^{J}$ with known realization.

2.1 Two-stage stochastic commodity production and supply planning

The application of commodity production and supply in an oligopolistic market serves as a motivation as well as the practical problem of interest. Presented as a stochastic game, the strategy of each agent in supply-side of the market can be described as the solution of a stochastic optimization problem (3). The decision process follows that agent $j\in\mathcal{J}$ decides an optimal production quantity $x_{j}$ of the commodity at the production stage. At the second-stage, each agent decides a supply quantity $y_{j}(\xi)$ after $\xi$ is observed, and a total quantity $T(y(\xi)):=\sum_{j=1}^{J}{y_{j}(\xi)}$ is supplied to the market. Our focus on oligopolistic markets requires that the price is dominantly affected by the total supplied quantity in the market $T(y(\xi))$. Therefore, all the trading occurs at the price $p:\Xi\times\mathbb{R}\to\mathbb{R}_{+}$, determined by a stochastic inverse demand curve $p(\xi,T(y(\xi))$. In practice, production and supply quantities are subject to physical restrictions, e.g., capability of production plant, logistic restriction, etc., i.e., $x_{j}\in X_{j}$ and ${y_{j}(\xi)}\in_{a.s.}Y_{j}$. More specifically, we have non-negative requirements for both production and supply, $X_{j}=\mathbb{R}_{+}$ and $Y_{j}=\mathbb{R}_{+}$. The relations between stage-wise decision variables $x_{j}$ and $y_{j}(\xi)$ are captured by constraints g_{j}\big{(}\xi,x_{j},{y_{j}(\xi)}\big{)}\leq_{a.s.}0 in (3). Essentially, every agent needs to formulate and solve a two-stage stochastic programming problem with recourse in the sense of achieving equilibria of a $J$-agents non-cooperative game of the market. We further require that agent $j$’s supply to the market cannot exceed his/her production quantity, i.e., ${y_{j}(\xi)}-x_{j}\leq_{a.s.}0$. This can be interpreted as the fact that agents may have no stock to start with, or they need to preserve certain reserved quantities prior to each decision process.

The problem can then be viewed from a slight different perspective than that of problem (3). As seen from the first-stage, agent $j\in\mathcal{J}$ wants to find a production quantity $x_{j}\geq 0$ to

[TABLE]

where,

[TABLE]

Objective function in (4) is regarded as the expected profit of agent $j$’s and problems (4)-(5) are termed intrinsic first-stage problem following that of a related treatment in [35]. In particular, the analysis of intrinsic first-stage problem and the stochastic programming problem with recourse in convex case were carried out in a series of studies by Rockafellar and Wets [34, 35, 36, 37] and more recently in [40]. The key feature of intrinsic first-stage problem, as well as formulation (3), is the requirement on precise orders of decision execution, commonly known as the constraints of nonanticipativity, which implies the decision variable in each stage should only depend on the information (data) available up to that stage. In (4)-(5), the second-stage decisions are explicitly determined after the first-stage decision, provided for each $x_{j}$ the second-stage problem is well-defined [43]. However, the study of optimality condition of (4)-(5), in the case of a general probability space $(\Xi,\mathcal{F},P)$333In cases of finitely supported distribution, the equivalence between intrinsic first-stage problem and the original recourse problem can be established, and the optimality condition of the recourse problem can be applied, see for example [40]., is very complicated since one needs to characterize the order of the decision process explicitly. For ease of analysis, we assume that there exists a multiplier $\lambda_{j}\in\mathfrak{L}^{1}_{1}$ corresponds to second-stage constraint and study the saddle-point condition of the Lagrangian formulation of problem (3). It is worth mentioning that the Karush-Kuhn-Tucker (KKT) condition of problem (3) (see [34]) introduces a second-stage multiplier $\tilde{\lambda}_{j}\in(\mathfrak{L}^{\infty}_{1})^{*}$ for every $j\in\mathcal{J}$ which incorporates the two-stage decision making process. This can be seen from the fact that any element of the dual space $(\mathfrak{L}^{\infty}_{1})^{*}$ can be decomposed into a component of $\mathfrak{L}^{1}_{1}$ and a “singular” component, corresponding to the multiplier of nonanticipativity. The saddle-point condition is shown to be sufficient and “almost” necessary for optimality of problem (3), and we refer the interested readers to [34, 35, 36, 37, 40] for more details.

The Lagrangian formulation of problem (3) associated with agent $j$ is of the following form

[TABLE]

where

[TABLE]

The constraints ${y_{j}(\xi)}\leq_{a.s.}x_{j}$ can be interpreted as the situation under which the profit maximizing supply $y^{*}_{j}(\xi)$ of agent $j$ is not necessarily equal to the total production quantity $x_{j}$. This feature of our model differs from conventional requirement on production-clearing condition, i.e., all the produced goods are expected to supply to the market.

In order to make further progress in characterizing the CN equilibrium, we need to specify the structures of our wealth functions, $W^{1}_{j}$ and $W^{2}_{j}$, suitable for our application. We assume that the production cost for $j$-th agent is quadratic, i.e., for each $j\in\mathcal{J}$ the cost of producing $x_{j}$ amount of production is $\frac{1}{2}c_{j}x^{2}_{j}+a_{j}x_{j}$, for some $c_{j}>0,a_{j}>0$. In the second-stage, the cost function of the supply or second-stage is assumed to be linear and of stochastic nature, i.e., for each $j\in\mathcal{J}$ the cost of supplying ${y_{j}(\xi)}$ amount of commodity is ${h_{j}(\xi)}{y_{j}(\xi)}$ for a.e. $\xi\in\Xi$. Here, $h_{j}(\xi)$ can be regarded as the unite supply cost of agent $j$. We adopt a classic stochastic inverse demand curve, see for example [23], that takes the expression p\big{(}\xi,T(y(\xi))\big{)}={p_{0}}(\xi)-\gamma(\xi)T(y(\xi)) for the spot price. In practice, the stochastic benchmark price excluding the effect of supply to the market $p_{0}:\Xi\mapsto\mathbb{R}_{+}$ can be estimated via statistical approaches based on real data. The supply discount $\gamma:\Xi\rightarrow\mathbb{R}_{+}$ acts as a market mechanism to adjust and reflect uncertainty in quantity in the market. In order to respect the market mechanism of supply-demand relation, we make the following assumption through out our study. {assumption} There exists a $\gamma_{0}>0$ such that $\gamma(\xi)\geq\gamma_{0}$ for a.e. $\xi\in\Xi$.

Thus, agent $j$’s stage-wise wealth functions are,

[TABLE]

and

[TABLE]

where the short-handed notation of the risk-adjusted spot price of agent $j$’s is denoted by $p_{j}(\xi):=p_{0}(\xi)-h_{j}(\xi).$

We are now ready to consider the specific stochastic programming for every agent $j\in\mathcal{J}$ :

[TABLE]

where

[TABLE]

Note that the requirements in problem (7) hold a.s. in accordance with the a.s. constraints of the second-stage in problem (3). However, (6)-(7) is not easy to solve, especially in a SNEP with $J\geq 2$, see [8]. The complication arises since the $j$-th agent’s problem contains that of the other agents’ strategy, not yet known at the decision horizon. A commonly used method is to recast problem (6)-(7) of each agent as a stochastic equilibrium problem. Then, obtaining an equilibrium of the convex $J$-player game (6)-(7) is equivalent to finding its solution for all agents. Stochastic equilibrium has been shown to be an effective method to study and to solve two-stage multi-players stochastic game problems, see for instance [9, 28, 46, 40, 33]. We study the saddle-point condition of the problem (6)-(7), rewritten in the form of problem (3). More specifically, for all $j\in\mathcal{J}$, there exists $\bar{\lambda}_{j}(\xi)\in\mathfrak{L}^{1}_{1}$ with ${\bar{\lambda}(\xi)}\geq_{a.s.}0$ so that a strategy \big{(}\bar{x}_{j},\bar{y}_{j}\big{)}\in\mathbb{R}_{+}\times\mathfrak{L}^{\infty}_{+} solves the following system.

[TABLE]

In particular, stationarity comes from the first-order necessary optimality condition under the assertion \partial\mathbb{E}\big{[}\Phi_{j}(\xi,x)\big{]}\subseteq\mathbb{E}\big{[}\partial_{x}\Phi(\xi,x)\big{]}. The assertion is discussed in [2], and the above system can be viewed as a weaker condition for optimality.

Rewritten in a compact form as SVI, the optimal strategy-multiplier pair $(x_{j},y_{j},\lambda_{j})\in\mathbb{R}_{+}\times\mathfrak{L}^{\infty}_{+}\times\mathfrak{L}^{1}_{+}$ must satisfy,

[TABLE]

It follows that since all agents in oligopolistic market act non-cooperatively, we write down the equilibrium interpreted as that of the whole system. More specifically, let $x=(x_{1},\ldots,x_{J})^{T}$ be the first-stage decision vectors of the system, and for a.e. $\xi\in\Xi$, the second-stage decision vector y(\xi)=\big{(}y_{1}(\xi),\ldots,y_{J}(\xi)\big{)}^{T} and its corresponding multiplier vector \lambda(\xi)=\big{(}\lambda_{1}(\xi),\ldots,\lambda_{J}(\xi)\big{)}^{T} are denoted respectively. Similarly, parameter vectors can be written collectively as $a=(a_{1},\ldots,a_{J})^{T}$, $p(\xi)=(p_{1}(\xi),\ldots,p_{J}(\xi))^{T}$. Then, we can treat the SVI for all agents as a two-stage stochastic complementary problem (SCP):

[TABLE]

where

[TABLE]

It follows that for the whole system, a $J$-tuple of strategies

[TABLE]

is called a solution of the two-stage SCP (9).

3 Structure of the regularized two-stage SCP

In this section, we focus on characterizing solutions of two-stage stochastic linear complementarity problem (9). From the derivation of first-order necessary optimality conditions of problem (6)-(7) and the monotonicity of problem (9), we have the following results on existence of solutions.

Proposition 3.1** (Theorem 2, [41]).**

For any fixed pair $(x,\xi)\in\mathbb{R}_{+}^{J}\times\Xi$, the second-stage problem (7) has a unique solution.

Thus, for the two-stage stochastic linear complementarity problem (9), the following proposition holds.

Proposition 3.2**.**

The two-stage stochastic linear complementarity problem (9) has relatively complete recourse, i.e., for any $x\in\mathbb{R}_{+}^{J}$ and a.e. $\xi\in\Xi$ the second-stage problem of (9) is solvable.

Proof 3.3**.**

The coefficient matrix of the second-stage part of (9)

[TABLE]

is positive semidefinite for a.e. $\xi\in\Xi$.

For any given $x\in\mathbb{R}_{+}^{J}$, it follows that there always exists a pair \big{(}\hat{y}(\xi),\hat{\lambda}(\xi)\big{)}\in\mathbb{R}^{J}\times\mathbb{R}^{J}, such that

[TABLE]

In detail, we consider a special choice $\hat{y}(\xi)=0$ and $\hat{\lambda}(\xi)=\max\{0,p(\xi)\}$, where the $\max$ function is taken componentwise. Thus, the corresponding quadratic programming problem of the linear complementarity problem is feasible. It follows from [32, Lemma 3.1.1, Theorem 3.1.2] that there must exist at least a solution which solves the second-stage problem for any given pair $(x,\xi)$. ∎

Although the second-stage problem (7) has a unique equilibrium for any given $(x,\xi)$ (see Proposition 3.1), the system (9) may admit multiple solutions. To see this, we give an illustrative example.

Example 3.4**.**

Consider a duopoly game, with given $x=(x_{1},x_{2})^{T}\geq 0$, and $-p(\xi)\geq_{a.s.}0$. Then, the corresponding second-stage part of complementarity system (9) reads

[TABLE]

Then, the solution set of (10) is of the following form

[TABLE]

where $\lambda_{1}(\xi)\geq_{a.e.}0,\lambda_{2}(\xi)\geq_{a.e.}0$.

In Example 3.4, the “equilibrium price” $\lambda$ may admit multiple values when there exist some zero-valued components of $x$.

Technically, the multiple solutions of the second-stage problem will cause trouble when we handle the two-stage stochastic complementarity system (9), both in computation and analysis [43]. The assumption ensuring the uniqueness of second-stage solution is usually made, see for instance [8, 9]. Moreover, interpreted as “equilibrium price” associated with agents’ production clearing, different values of $\lambda$ would have ambiguous economical interpretations. Motivated by these, we propose a regularized method to seek for one particular choice of “equilibrium price”. Similar approach can be found in for example [7].

For an $\epsilon>0$, let

[TABLE]

Thus, we propose the regularized SCP of (9) as follows:

[TABLE]

For a given pair $(x,\xi)\in\mathbb{R}_{+}^{J}\times\Xi$, the second-stage problem of (9) and the regularized second-stage problem (11) are denoted by LCP$(q(x,\xi),M(\xi))$ and LCP$(q(x,\xi),M^{\epsilon}(\xi))$ respectively. Their solution functions are chosen from the respective solution sets and expressed by $z(q(x,\xi))$ and $z^{\epsilon}(q(x,\xi))$. In the sequel, we omit the $\xi$ and $x$ without causing confusion, i.e., LCP$(q,M):=\text{LCP}(q(x,\xi),M(\xi))$ and LCP$(q,M^{\epsilon}):=\text{LCP}(q(x,\xi),M^{\epsilon}(\xi))$.

For clearer demonstration, recall our illustrative Example 3.4, and consider its regularization approach. Thus, the second-stage of the regularized problem takes the following form

[TABLE]

Under the same condition as in Example 3.4, we can obtain the unique solution of (12), which $\tilde{y}_{1},\tilde{y}_{2},\tilde{\lambda}_{1},\tilde{\lambda}_{2}$ equal to 0 for a.e. $\xi\in\Xi$. Due to the positive definiteness of $C$, it follows that we obtain the unique solution of the first-stage problem is $x_{1}=0,x_{2}=0$. Then, we have obtained one particular solution of the original problem, the trivial solution in this example. The key feature of our regularized method is that it promises the existence and uniqueness of solution due to the strongly monotone of regularized two-stage problem.

In the remaining of this section, we concern ourselves with the solution $z^{\epsilon}$ of LCP$(q,M^{\epsilon})$ and explore the structure of the second-stage solution.

Proposition 3.5**.**

For any fixed $\epsilon>0$, the regularized problem (11) has a unique solution $(x^{\epsilon},y^{\epsilon},\lambda^{\epsilon})\in\mathbb{R}^{J}\times\mathcal{Y}\times\mathcal{Y}$.

Proof 3.6**.**

The result can be obtained via a similar procedure as in [9, Proposition 2.1 (i)] and we only need to show that the condition in [9, Assumption 1] holds. Recall that Assumption 2.1 holds, then for a.e. $\xi\in\Xi$

[TABLE]

where $\tau=\min\{\bar{c},\gamma_{0}(J+1),\epsilon\}$ with $\bar{c}$ denoting the minimum diagonal element of $C$. ∎

Theorem 3.7**.**

For any fixed $\epsilon>0$, $x\geq 0$ and a.e. $\xi\in\Xi$, the $j$-th component of the second-stage solution of problem (11) \big{(}(y^{\epsilon})_{j},(\lambda^{\epsilon})_{j}\big{)} is either $(0,0)$, or one of the following two forms:

[TABLE]

for $j\in\mathcal{J}$, where

[TABLE]

with

[TABLE]

where $|\mathcal{I}_{2}|$ and $|\mathcal{I}_{3}|$ denote the cardinality of $\mathcal{I}_{2}$ and $\mathcal{I}_{3}$ respectively.

The detailed proof of the Theorem 3.7 is given in Appendix. Note that the above theorem gives the forms of the unique solution of the second-stage regularized problem (11). However, it may be used to assist numerical calculation since the partition of the index set is not known in advance. Nevertheless, it is suffice for our purposes of deriving additional properties of the solutions. Due to the positive definiteness of $M^{\epsilon}$ and special structure of problem (11), we first obtain the following Lipschitz continuous property, following [11, Corollary 2.1].

Lemma 3.8**.**

For any $\xi\in\Xi$, there exists $L(\xi)>0$ such that for any fixed $\epsilon\in(0,1]$, we have

[TABLE]

Lemma 3.9**.**

For any fixed $\epsilon>0$ and $(x,\xi)\in\mathbb{R}_{+}^{J}\times\Xi$, $T^{\epsilon}$ has the following upper bound:

[TABLE]

Proof 3.10**.**

We have the following derivation from (14) that

[TABLE]

∎

We end this section by establishing the convergence result of the second-stage LCP$(q,M^{\epsilon})$ solutions as $\epsilon\downarrow 0$ for any given pair $(x,\xi)\in\mathbb{R}^{J}\times\Xi$.

Proposition 3.11**.**

For any fixed $\epsilon>0$ and $(x,\xi)\in\mathbb{R}_{+}^{J}\times\Xi$, let $z^{\epsilon}(\xi)=(y^{\epsilon}(\xi),\lambda^{\epsilon}(\xi))$ denote the unique solution of the regularized problem LCP$(q,M^{\epsilon})$. Then

[TABLE]

where $\bar{z}(\xi)=(\bar{y}(\xi),\bar{\lambda}(\xi))$ denotes the unique least $l_{2}$-norm solution of the LCP$(q,M)$. Moreover, the $j$-th component of the least $l_{2}$-norm solution of problem (9) has one of the following three forms:

[TABLE]

for $j\in\mathcal{J}$, where

[TABLE]

Furthermore, for a.e. $\xi\in\Xi$ there exists $\bar{\kappa}(\xi)>0$, such that

[TABLE]

The proof of Proposition 3.11 is given in Appendix.

4 Convergence analysis

In this section, we first prove the convergence of the unique solution of the regularized problem (11) to the solution set of the original problem as the regularized parameter $\epsilon$ decreases to zero. Next, we will study the sample average approximation (SAA) to solve the regularized problem. See [7]. Combined with our regularization approaches, we demonstrate the convergence property of the solution of our regularized SAA model as the number of samples goes to infinity. More specifically, the convergence analysis in this section is divided into two parts: the convergence analysis of the regularized problem as the regularized parameter $\epsilon$ tends to zero, and the analysis of regularized SAA. We finally build up the convergence relationship between the regularized SAA approach and the original problem.

4.1 Convergence of the regularized model

In this subsection, we only need to consider the convergence properties of the first-stage decision vector, i.e., $x^{\epsilon}\in\mathbb{R}^{J}_{+}$ that solves problem (11), when the regularized parameter $\epsilon$ tends to zero. The convergence property of the solution $(x^{\epsilon},y^{\epsilon},\lambda^{\epsilon})$ then follows combining the result of section 3. From Proposition 3.5, we know that for fixed $\epsilon>0$ problem (11) admits a unique first-stage solution $x^{\epsilon}$. In the following, we concern about the sequence of accumulation points of $\{x^{\epsilon}\}$ as $\epsilon\downarrow 0$.

For the existence of accumulation points, we have the following result.

Proposition 4.1**.**

Suppose there exists $p_{0}>0$ such that for all $j\in\mathcal{J}$, $p_{j}(\xi)\leq_{a.s.}p_{0}$. Then, with $\epsilon\downarrow 0$, $\{x^{\epsilon}\}$ is bounded.

Proof 4.2**.**

From the condition on $p_{j}(\xi)$, there must a sufficiently large $\alpha>0$ such that for any $j\in\mathcal{J}$

[TABLE]

Then we have

[TABLE]

Assume that $\{x^{\epsilon}\}$ is unbounded for the purpose of arriving at a contradiction. Then, it follows that there exist some indices $j\in\mathcal{J}$, such that $(x^{\epsilon})_{j}\geq\alpha$. Then we consider the $j$-th component of the first-stage complementarity relation,

[TABLE]

which can be expressed, from (13), as

[TABLE]

However, this complementarity relation cannot be obtained because $(x^{\epsilon})_{j}>0$ and $c_{j}(x^{\epsilon})_{j}+a_{j}>0$. This completes our proof. ∎

Note that the conditions $p_{j}(\xi)\leq_{a.s.}p_{0}$ can be easily satisfied in many practical applications. For example, with given data sets of $p(\xi)$ we can always find an upper bound $p_{0}:=\max_{j}\{p_{j}(\xi)\}$.

Lemma 4.3**.**

Suppose there exists a constant $p_{0}>0$ such that for all $j\in\mathcal{J}$, $p_{j}(\xi)\leq_{a.s.}p_{0}$. Then, there exists a sequence $\{\epsilon_{k}\}_{k=1}^{\infty}$ with $\epsilon_{k}\downarrow 0$ as $k\rightarrow\infty$ such that $x^{k}\to\hat{x}$ and $\lambda^{k}\to\hat{\lambda}$,

[TABLE]

where $x^{k}$ and $\lambda^{k}(\xi)$ are parts of the unique solution of problem (11) with $\epsilon=\epsilon_{k}$, and $\bar{\lambda}(\xi)$ is part of the least norm solution of the second-stage problem (11) with $x=\hat{x}$.

Proof 4.4**.**

From Proposition 4.1, there exists a sequence $\{x^{k}\}$ such that $\lim_{k\to\infty}x^{k}=\hat{x}$. From (16), we have

[TABLE]

From Lemma 3.8 and (16), for a.e. $\xi\in\Xi$

[TABLE]

We have from the derivation that the estimation

[TABLE]

where the last term comes from Lemma 3.8 with some vector $w$ with $\{x^{k}\}_{k=1}^{\infty},\hat{x}\subset[0,w]$. It follows from the Lebesgue Dominated Convergence Theorem, we have

[TABLE]

Then we complete the proof. ∎

Theorem 4.5**.**

Any accumulation point of $\{x^{\epsilon},y^{\epsilon},\lambda^{\epsilon}\}$ as $\epsilon\downarrow 0$ is a solution of problem (9).

Proof 4.6**.**

We only need to verify that for any $\epsilon_{k}\downarrow 0$, the accumulation point $\hat{x}$ of subsequence $\{x^{\epsilon_{k}}\}$ is a first-stage solution of (9). Since $x^{\epsilon_{k}}$ is the first-stage solution of problem (11) for any $\epsilon_{k}>0$, we have with $x^{k}=x^{\epsilon_{k}}$

[TABLE]

which, by using the ‘min’ NCP function (see, for example, [32]), can be rewritten as

[TABLE]

By Lemma 4.3, we have

[TABLE]

as $k\rightarrow\infty$. Thus we obtain that

[TABLE]

The statement then follows from Proposition 3.11. ∎

4.2 Convergence of the regularized SAA model

In this subsection, we study the SAA scheme for solving the regularized problem (11) and focus on the convergence of the regularized SAA approach. More specifically, we focus on the SAA convergence analysis and the solution of the first-stage problem. It is noteworthy that Chen, Sun and Xu considered a discrete approximation scheme in [9], which also leads to an approximation of the response variable in the second-stage problem.

Let $\xi_{1},\xi_{2},\ldots,\xi_{\nu}$ denote $\nu$ independent identically distributed (i.i.d.) samples. Then, with slight abuse of notation, we can obtain the following formulation of problem (11) with SAA:

[TABLE]

Or, we can write the problem collectively for all $\nu$ samples

[TABLE]

where $C\in R^{J\times J},$ B=\left(\begin{array}[]{cc}0&-I\\ \end{array}\right)\in R^{J\times 2J}, and for $\ell=1,\ldots,\nu$, D^{\epsilon}_{\ell}=\left(\begin{array}[]{cccc}\Pi(\xi_{\ell})&I\\ -I&\epsilon I\\ \end{array}\right)\in R^{2J\times 2J}, $v_{\ell}=(y(\xi_{\ell}),$ $\lambda(\xi_{\ell}))^{T},q_{\ell}=(-p(\xi_{\ell}),0)^{T}.$ Thus, (18) is treated as a large-scale deterministic linear complementarity problem:

[TABLE]

where $z=(x,v_{1},...,v_{\nu})^{T}$, $\bar{q}=(a,q_{1},\cdots,q_{\nu})^{T}$, and $H^{\epsilon}$ denotes the coefficient matrix in (18).

We have the following assertion of existence and uniqueness of problem (17) by [32, Theorem 3.1.6].

Proposition 4.7**.**

For any fixed $\epsilon>0$ and positive integer $\nu$, there exists a unique solution of problem (17).

Recall the result of Lemma 3.8 and the following proposition can be shown in a similar way as in [7, Proposition 3.7].

Proposition 4.8**.**

Let $(y^{\epsilon}(\xi),\lambda^{\epsilon}(\xi))$ be the unique solution of the regularized second-stage problem (11) for any $(x,\xi)\in\mathbb{R}_{+}^{J}\times\Xi$. Then,

[TABLE]

with probability (w.p.) $1$ as $\nu\rightarrow\infty$ uniformly on $\mathbb{B}(x,\delta)\cap\mathbb{R}_{+}^{J}$ for any $\delta>0$,

Let $x^{\epsilon}_{\nu}$ denote the first J-components of the unique solution of problem of (17), and we have the following assertion.

Lemma 4.9**.**

Suppose there exists $p_{0}>0$ such that for all $j\in\mathcal{J}$, $p_{j}(\xi)\leq_{a.s.}p_{0}$. Then, with $\epsilon\downarrow 0$, $\{x_{\nu}^{\epsilon}\}$ is bounded.

We omit the proof since it can be shown analogously as in Proposition 4.1.

Theorem 4.10**.**

Suppose there exists $p_{0}>0$ such that for all $j\in\mathcal{J}$, $p_{j}(\xi)\leq_{a.s.}p_{0}$. Then, for any fixed $\epsilon>0$, $x^{\epsilon}_{\nu}\rightarrow x^{\epsilon}$ w.p. $1$ as $\nu\rightarrow\infty$.

Proof 4.11**.**

From Propositions 3.5, Proposition 4.7, and Lemma 4.9, for any fixed $\epsilon>0$, both the regularized problem (11) and its SAA-regularized problem (17) have solutions and contained in some compact subset in $\mathbb{R}_{+}^{J}$. We know from Proposition 4.8 that

[TABLE]

as $\nu\rightarrow\infty$, uniformly with respect to $x$ on any compact set. Then, we have $x^{\epsilon}_{\nu}\rightarrow x^{\epsilon}$ w.p. $1$ as $\nu\rightarrow\infty$ by [44, Proposition 19]. ∎

Combining Theorem 4.5 with Theorem 4.10, we have the following convergence result.

Theorem 4.12**.**

Suppose there exists $p_{0}>0$ such that for all $j\in\mathcal{J}$, $p_{j}(\xi)\leq_{a.s.}p_{0}$. Then,

[TABLE]

w.p. $1$, where $S^{*}$ denotes the optimal solution set of the first-stage problem of (9).

5 Numerical tests and applications on crude oil market

In this section, we firstly carry out numerical experiments using randomly generated data to illustrate the effectiveness of our model and its solution approach. Furthermore, we adopt our two-stage stochastic CN equilibrium problem to describe market share competition in crude oil market. The results show that the model is capable of reproducing the actual oil market share based on in-sample data. Moreover, it is also shown that the model can make good out-of-sample predictions using historical data. All the tests are run in MATLAB 2016b on a personal computer with 32GB RAM and 8-core processor ($3.6\times 8GHz$).

5.1 Progressive hedging method and smoothing Newton sub-algorithm

Firstly, randomly generated problems are used for testing our regularized SAA approach to solve the two-stage stochastic CN equilibrium problem. Recall that the model of interests takes the form of a scenario-based linear complementarity problem (17) or its equivalent expression (18) with sufficiently small $\epsilon$. The solution process adopts the well-known PHM. The PHM is globally convergent and the convergence rate is linear for problem (18) with SAA approach.

The classical PHM to solve (18) is as follows

Algorithm 1: Progressive hedging method

Step 0. Given an initial point $x^{0}\in\mathbb{R}^{J}$, let $x_{\ell}^{0}=x^{0}\in\mathbb{R}^{J},v_{\ell}^{0}\in\mathbb{R}^{2J}$ and $w_{\ell}^{0}\in\mathbb{R}^{J}$, for $\ell=1,\ldots,\nu$, such that $\frac{1}{\nu}\Sigma_{\ell=1}^{\nu}w_{\ell}^{0}=0.$ Set the initial point $z^{0}=(x^{0},v_{1}^{0},\ldots,v_{\nu}^{0})^{T}$. Choose a step size $r>0.$ Set $k=0.$

Step 1. If the point $z^{k}$ satisfies the condition

$\|\min(z^{k},H^{\epsilon}z^{k}+\bar{q})\|\leq 10^{-6},$

output the solution $z^{k}$ and terminate the algorithm; otherwise, go to Step 2.

Step 2. For $\ell=1,\ldots,\nu$, find $(\hat{x}^{k}_{\ell},\hat{v}_{\ell}^{k})$ that solves linear complementarity problems

\begin{array}[]{l}0\leq x_{\ell}\bot Cx_{\ell}+Bv_{\ell}+a+w_{\ell}^{k}+r(x_{\ell}-x_{\ell}^{k})\geq 0,\\ 0\leq v_{\ell}\bot-B^{T}x_{\ell}+D^{\epsilon}_{\ell}v_{\ell}+q_{\ell}+r(v_{\ell}-v_{\ell}^{k})\geq 0.\end{array}

(20)

Then let $\bar{x}^{k+1}=\frac{1}{\nu}\sum_{\ell=1}^{\nu}\hat{x}_{\ell}^{k}$, and for $\ell=1,\ldots,\nu,$ update

$x_{\ell}^{k+1}=\bar{x}^{k+1},v_{\ell}^{k+1}=\hat{v}_{\ell}^{k},w_{\ell}^{k+1}=w_{\ell}^{k}+r(\hat{x}_{\ell}^{k}-x_{\ell}^{k+1}),$

to get point $z^{k+1}=(\bar{x}^{k+1},v_{1}^{k+1},\ldots,v_{\nu}^{k+1})^{T}$.

Step 3. Set $k:=k+1;$ go back to Step 1.

The PHM involves solving $\nu$ independently sample-based LCP (20) at each iteration. Problem (20) is well-defined, since for $\ell=1,\ldots,\nu,$ the coefficient matrix

[TABLE]

is positive definite for any $\epsilon>0$. Thus, it has a unique solution for each $\ell=1,\ldots,\nu$. For simplicity, denote (20) as

[TABLE]

where z_{\ell}=\left(\begin{array}[]{c}x_{\ell}\\ v_{\ell}\\ \end{array}\right),\tilde{H}_{\ell}^{\epsilon}=\left(\begin{array}[]{cc}C+rI&B\\ -B^{T}&D^{\epsilon}_{\ell}+rI\\ \end{array}\right),\tilde{q}_{\ell}=\left(\begin{array}[]{cc}a+w_{\ell}^{k}-rx_{\ell}^{k}\\ q_{\ell}-rv_{\ell}^{k}\end{array}\right) with $\ell=1,2,\ldots,\nu.$ Then, we can equivalently write a large-scale LCP for all $\nu$ samples and solve

[TABLE]

The structure of problem (22) enables us to use block computation to solve it, which can significantly improves the efficiency of the PHM. For example, suppose $\nu=mN$ with $m,N$ being positive integer number. The equivalent $m$-block reformulation for (22) reads

[TABLE]

where $i_{j}=iN+j,j=1,\ldots,m,i=0,\ldots,(m-1)N$.

The main computation cost of the PHM is in Step 2 due to the large sample size $\nu$ despite the relative low cost in solving for one sample. Block implementation speeds up the computation by better exploring and rebalancing the computational load. In practice, the number of blocks are adjusted so that the over all computation time is minimized. To improve the efficiency of the PHM, we also use the warm-start technique suggested in [33] to choose an initial point for subproblem (20). More specifically, the solution $z^{k}$ of subproblem (20) at the $k$-th iteration is used as a starting point for the $(k+1)$-th iteration.

In the remaining of this subsection, we focus on problem (21) for a given sample and omit the subscription ℓ for ease of expression. In order to take advantage of the sparse structure of the subproblem (21), we apply the smoothing Newton method proposed by Chen and Ye [12] to solve it. In what follows, we give a brief description of the smoothing Newton method. It is well-known that solving (21) for a given sample is equivalent to solving the nonsmooth equation

[TABLE]

The main idea of the smoothing Newton method is to use a smooth approximation function to approximate the nonsmooth function $F$ and then solve the corresponding linear system. We use the smooth Gariel-Moré approximation function $f:\mathbb{R}^{3J}\times\mathbb{R}_{++}\rightarrow\mathbb{R}^{3J}$ to approximate the nonsmooth function $F$. The density function is $\rho(s)=2/(s^{2}+4)^{\frac{3}{2}}$ (see [12] for details). In our numerical tests, the $j$-th component of $f$ reads

[TABLE]

where the corresponding $j$-th diagonal element of the Jacobian $\bar{D}(z)$ is

[TABLE]

Then, the smoothing Newton method for solving subproblem (24) requires to solve a linear equation to determine $d^{k}$ at each iteration, namely

[TABLE]

where $\nabla_{z}f(z^{k},\delta_{k})=I-\bar{D}(z^{k})(I-\tilde{H}^{\epsilon})$, $\delta_{k}$ decreases to 0 according to the criterion in article [12]. To guarantee the well-defineness of (25), we make use of the following result.

Theorem 5.1**.**

[15]** For any diagonal matrix $\tilde{D}=diag(\tilde{D}_{jj})\in\mathbb{R}^{J\times J}$ with $0\leq\tilde{D}_{jj}\leq 1,j=1,2,\ldots,J$, the matrix $I-\tilde{D}(I-A)$ is nonsingular if and only if $A$ is a P-matrix. ∎

For any given sample, it is known that $H^{\epsilon}$ is positive definite and hence a P-matrix. Moreover, $\bar{D}(z)$ is a diagonal matrix with its element on the interval $[0,1]$ for any $(z,\delta)\in\mathbb{R}^{3J}\times\mathbb{R}_{++}$. Therefore, using Theorem 5.1, the Jacobian $\nabla_{z}f(z,\delta)$ is nonsingular for any $(z,\delta)\in\mathbb{R}^{3J}\times\mathbb{R}_{++}$. Thus, the linear equation (25) is well-defined.

Denoting the matrix $\bar{D}(z)=diag(\bar{D}_{1}(z),\bar{D}_{2}(z),\bar{D}_{3}(z))$, the Jacobian $\nabla_{z}f(z,\delta)$ at the point $z$ is of the following structure

[TABLE]

where $\Lambda_{1}(z)=\bar{D}_{1}(z)(C+(r-1)I)+I$, $\Lambda_{2}(z)=(\gamma(\xi)+(r-1))\bar{D}_{2}(z)+I$, $\Lambda_{3}(z)=(\epsilon+(r-1))\bar{D}_{3}(z)+I$ are all diagonal matrices, and $u_{1}(z)=\gamma(\xi)\bar{D}_{2}(z)e$. We can take advantage of the sparse structure of (29) in solving (25). Specifically, $\nabla_{z}f(z,\delta)$ consists of only diagonal sub-matrix and the matrix $u_{1}(z)e^{T}+\Lambda_{2}(z)$, where the later is a sum of a diagonal sub-matrix and a rank-one matrix.

Noticing from (29), linear equation (25) is of the following form

[TABLE]

where $\Lambda_{i}\in\mathbb{R}^{J},i=1,2,\ldots,7$ are diagonal matrices with $\Lambda_{1}$ and $\Lambda_{3}$ being nonsingular, $s_{i}\in\mathbb{R}^{J},b_{i}\in\mathbb{R}^{J},i=1,2,3,$ and $u_{1},u_{2}\in\mathbb{R}^{J}.$ For ease of notation, we use $\hat{\Lambda}=diag(1/\Lambda_{11},\ldots,1/\Lambda_{JJ})$ to represent the inverse of any invertible diagonal matrix $\Lambda=diag(\Lambda_{11},\ldots,\Lambda_{JJ})$. Given the sparse structure of the coefficient matrix, we can solve (30) efficiently. More exactly, by the first two equations of (30), we can get

[TABLE]

Directly substituting (31) and (32) into the third equation of (30), then we have

[TABLE]

where $const=-(\Lambda_{5}\hat{\Lambda}_{1}b_{1}+\Lambda_{6}(u_{1}u_{2}^{T}+\Lambda_{3})^{-1}b_{2})$. For computing the inverse matrix of $(u_{1}u_{2}^{T}+\Lambda_{3})$, the Sherman-Morrison formula is useful. We have by the Sherman-Morrison formula

[TABLE]

Substituting (34) into (33), we get

[TABLE]

where $\alpha=1/(1+u_{2}^{T}\hat{\Lambda}_{3}u_{1}),\Lambda_{0}=\Lambda_{7}-\Lambda_{5}\hat{\Lambda}_{1}\Lambda_{2}-\Lambda_{6}\hat{\Lambda}_{3}\Lambda_{4},\tilde{u}_{1}=\Lambda_{6}\hat{\Lambda}_{3}u_{1},\tilde{u}_{2}=\Lambda_{4}\hat{\Lambda}_{3}u_{2}.$ Then, if $\Lambda_{0}$ is nonsingular, using the Sherman-Morrison formula again, we can immediately get the solution of $s_{3}$

[TABLE]

Then, substituting the $s_{3}$ into (31) and (32), we get the solution of $s_{1}$ and $s_{2}$, respectively.

From (35), one can know that the computation cost of $s_{3}$ is trivial, since we only need to compute inverse of several diagonal matrix, namely, $\Lambda_{0},\Lambda_{1},$ and $\Lambda_{3}$. Once $s_{3}$ is obtained, the calculation of $s_{1}$ and $s_{2}$ just needs to perform matrix-vector production. Therefore, the linear equation $\eqref{smooth-linear-eq}$ can be solve efficiently.

5.2 Randomly generated problems

For the first part of numerical test, we randomly generated the problem of the form (18). More specifically, we generate a set of i.i.d. samples $\{\xi_{\ell}\}_{\ell=1}^{\nu}$ from a uniformly distribution over the interval $[0,1]$. For $\ell=1,\ldots,\nu$, set

[TABLE]

Diagonal matrix $C\in\mathbb{R}^{J\times J}$ and $a\in\mathbb{R}^{J}$ are generated with its elements uniformly distributed over the interval $[1,2]$. All the numerical results are based on the average of 10 independent runs.

To show the feasibility of the solution of the regularized problem compared to that of the original problem, we compute the following residual value

[TABLE]

where $H$ denotes the coefficient matrix of (18) with $\epsilon=0$.

Selected numerical results for $J=10$ were listed in the Table 1 and 2. The first table is the solution obtained by PHM without block computation, while table 2 illustrates results with block implementation. For our generated experiments, we found $N=50$ is the best block choice (measured by CPU time) for $J=10$. The average number of iterations, the average cpu time, and the average value of Res were recorded in both table. It is easy to see that the block implementation greatly reduces the cpu time. For the same value of $\epsilon$, the number of iterations increases slightly when the sample size $\nu$ increases. In cases where the sample size $\nu$ is kept constant and the values of regularization parameter $\epsilon$ are chosen from $\epsilon=10^{-3}$ to $\epsilon=10^{-12}$, the iteration numbers are barely influenced as well as the cpu time. Furthermore, we observe the convergence of our regularization approach with decreasing values of $\epsilon$, as have been proved in previous sections. Also notice that, the value of Res decreases when the $\epsilon$ diminishes from $10^{-3}$ to $10^{-12}$. Numerically, it shows that the solution of the regularized problem is also that of the original problem when $\epsilon=10^{-12}$.

Figure 1 illustrates the performance of the PHM measured by the number of iterations and cpu time for $10$ players. It is also worth mentioning that although one might expect the problem to be more difficult to solve for a small $\epsilon$, the numerical performance in our experiments remain roughly unaffected with decreasing values of $\epsilon$.

Figure 2 demonstrates the convergence property of the first-stage solution $x$ when the sample size gets large for the case $J=10$. The convergence trend can be seen component-wisely as the solution $x$ converges when the sample size $\nu$ gets large.

5.3 Crude oil market

In this subsection, we demonstrate one application of our two-stage stochastic Nash equilibrium model in describing the crude oil market share. Namely, we investigate the strategies of crude oil exporting agents via solutions of our reformulated SVI problem so as to recreate actual market shares. Based on historical data on crude oil market, we make in-sample back-tracking test to establish the effectiveness and validity of our model while explaining the market behaviour. Furthermore, the out-of-sample prediction capability of our model is demonstrated when training data is used to specify model parameters. From the results of our numerical tests, we conclude that our model is suitable to reproduce, predict and potentially capable to explain stable market shares of crude oil.

Crude oil market is one of the most widely studied commodity market in the world. The market has long been described as of an ‘intermediate between monopoly and perfect competition [22]”, supported by extensive historical data and market observation. When the market encounters large sudden events, it took “oil shocks” [17, 18, 25]. Conceptually, consequences of supply and demand fluctuation, along with many other factors, are reflected most directly in dramatic changes in oil price. One interesting observation is that the market share behaves rather smoothly even during periods of oil shocks [25]. Majority of the world’s crude oil is supplied by a few large oil exporting countries and they are viewed collectively as a finite number of large agents from which price-taking consumers purchase product at the same price [45]. Hence, it is no surprise that CN approaches have been adopted from earlier days of study in this field, see [42]. The study of oil shock gives rise to rich field of economic researches in oil market structure and as well as factors affecting oil price, see [26, 49, 1], etc. Our model does not attempt to make explanation or future predictions on oil shocks but rather aims to make sense of the stable characteristics of supply-side of oil market share: major agents acting non-cooperatively to achieve market equilibrium. We focus on the supply factor of oil, while treat other factors, e.g., demand, world economic situation, population, etc., as known information with uncertainty. Our model is found to be able to reproduce the market shares as well as to forecast future production plans of different oil-producing countries based on historical data. In short, we firstly use historical data to determine model parameters as well as to approximate distribution for uncertainty in observations of the market. Model parameters are tested within a fixed sampling window so that the in-sample results obtained by solving our model matches that of training data set. It then follows that over a short decision horizon in the future, we can adopt the trained model to prediction future production plans. More specifically, we took i.i.d. samples using a fixed-size rolling window sampling method while the prediction is obtained by solving our model with the adjusted set of parameters obtained by in-sample tests within the training window.

5.3.1 Data set, parameter selection and uncertainty description

In order to implement real market data into our model (6)-(7), we need to determine parameters $c_{j},a_{j}$ correspond to production costs of $j$-th agent by analyzing training set. In practice, model parameters are adjusted via a brute force learning process. More specifically, if the training set consists only one trading day of information, we adjust the parameters so that the solution, i.e., production quantities of $J$ agents, matches the estimated historical observation. When more data become available, the parameters are adjusted so that the average of solutions is close to the real estimate. It is found from our numerical implementation that the best fit solutions are obtained when values of $c_{j}$ and $a_{j}$ are taken to be inversely proportional to known market share within the sampling period. It is firstly learned that the common belief of “lower unit cost would lead to large sale” does not hold in our model, which also motivates the choice of quadratic production cost in our model. Moreover, it is found in our numerical experiment that the smaller the quadratic cost parameters $c_{j}$ the greater market share the agent would have. This can be understood from the fact that the large oil producing agents are more flexible in adjusting their production quantities as compare to smaller agents. A more mathematical interpretation is to think of $c_{j}$ as related to the penalty parameters of the augmented Lagrangian formulation corresponding to production constraints of each agent: production cost is upper bounded with some unit production cost $\tilde{a}_{j}$. In this way, one can think of this observation been that the countries has greater reserve of oil suffers less when the production upper bound gets violated. Nevertheless, we found that the quadratic is more than adequate to reproduce historical observations.

The scenarios of stochastic parameters are taken from empirical distribution of the sampling data set, guided by the analysis of [19]. In detail, we draw samples to generate scenarios of market’s supply discounting factor $\gamma(\xi)$ and $j$-th agent’s adjusted price $p_{j}(\xi)=p_{0}(\xi)-h_{j}(\xi)$. The data set used in this study are from (i)Statistical Review of World Energy444https://www.bp.com/en/global/corporate/energy-economics/statistical-review-of-world-energy.html, yearly data from 1984, published every June by bp. Inc; (ii) Oil Price Dynamics Report555https://www.newyorkfed.org/research/policy/oil_price_dynamics_report, weekly by Federal Reserve Bank of New York since 1986; (iii)U.S. Energy Information Administration666https://www.eia.gov, weekly and daily spot price of Brent. The empirical distribution are generated from the contents of oil dynamic report, consists of time series of weekly percentage change of Brent spot price and its corresponding components’ contributions, i.e., those from demand, residual, and supply. More specifically, the percentage change in price in $(k+1)$-th777$k$ is used to denote the order in time, distinguished from that of i.i.d. samples week compared to that of the $k$-th week is denoted by $\triangle p^{k}$, while the contribution of demand and supply changes are denoted as $\triangle D^{k}$ and $\triangle S^{k}$ respectively. The contribution to the price change that does not correspond to demand and supply is denoted as residual $\triangle R^{k}$. Therefore, for any week within the report, it holds that the Brent price change is determinately represented by

[TABLE]

Recall that our expression for (scenario-based) inverse demand function of quantity offered to the market for a given scenario $\xi_{\ell}$ is in the following form

[TABLE]

In accordance with data structure, $p^{k}(\xi_{\ell})$ corresponds to the adjusted predicted price from contributions other than supply, i.e.,

[TABLE]

where $p_{0}^{k}$ denotes the known Brent spot price prior to that of the concerned time. Both $\triangle d^{k}(\xi_{\ell})$ and $\triangle r^{k}(\xi_{\ell})$ are random scenarios taken from empirical distributions of historical demand and residual distributions within the sampling set respectively. To estimate $\gamma^{k}(\xi)$, corresponding to $k$-week data, we need the data of spot quantity supplied to the market, $T^{k}$. In practice, it is very difficult to obtain reliable data on total supply to the market over short observation window. Rather, we observe that trust-worthy estimate on daily oil supply based on annual data is available along with the observation that a steady growth of about 1% per year over the last four decades, according to Statistical Review of World Energy. Therefore, we take random instances of daily market supply $T^{k}(y_{\xi})$ which are taken from a uniformly distributed interval between 99% and 101% of yearly based daily estimate. Then, we can generate a set of data of stochastic excessive supply discount factor $\{\gamma^{k}(\xi_{\ell})\}$,

[TABLE]

where $p_{0}^{k}$ is known with certainty within the testing data set, and the $|\cdot|$ ensures that increase in quantity has a negative influence on price.

We present the results for our numerical experiments on reproduced oil market share. As shown in Fig 4, the in-sample results over the periods of oil shocks 2007, 2009 and 2014 are reproduced using solely from the available information within that year. For example, the in-sample experiment of 2009 uses all the available data obtained within 2009, e.g., daily price from 01/01/2009 to 31/12/2009. For the generation of random instances in-sample, we use weekly price contribution covering the year 2009 to form empirical distributions for supply excluded spot prices $p(\xi_{k})$, and discounting factor $\gamma(\xi_{k})$, from which we took i.i.d. samples to represent random scenarios. Hence, for in-sample experiments, our goal is to reproduce known market share with the adjustments of cost parameters $a_{j},c_{j}$ corresponding to $j$-th oil producing agent.

The out-of-sample tests cover the period from 2007 to 2017 with the one year length equivalent of sampling window in Figure 4. Note that one year sampling window is adopted since independent numerical results of selected years demonstrated that the market has a short memory in the sense that longer historical data would not provide extra information, at least in the setup of our model and approaches. Prior to the date of interests, all the historical data are assumed to be available while the rolling window sampling permits the usage of newly acquired data as soon as it becomes available. From the construction and interpretation of our two-stage model, it is reasonable to assume that cost of oil production remain unchanged over a short decision horizon. Furthermore, we tested different lengths of the sampling window, especially during times of oil shocks for better in-sample fitting, and the best fit sampling window length is used in out-of-sample tests. For example, we are interested in obtaining the market share of 2009, during times of recent financial crisis, using the fixed length sampling windows whose length varies from one year to five years. If we want to predict the production on 01/01/2009 using a sampling window of one year, we took one year sampling window from 01/01/2008 to 31/01/2008 to obtain empirical distributions from which scenarios of $\gamma(\xi_{k})$, $p(\xi_{k})$ are sampled. Taking $\nu$ samples, the solutions method follows that described in previous subsection. For the result on 02/01/2019, the sampling window rolls forward using the data form 02/01/2008 to 01/01/2009 since the last day information becomes known, while the model remains unchanged. Therefore, for the market share prediction of any given year, we obtained the daily results by taking average of 250 samples, while the yearly estimate is formed by taking the average on daily results. From Figure 4, it is observed that we can obtain fair good prediction quantities through out the testing period. However, it is worth mentioning that during the numerical experiments across that decision horizon the computational cost is much higher during periods over oil shocks compare to those of relative stable oil prices.

6 Conclusion

A two-stage SVI is used to describe the equilibrium of a convex two-stage non-cooperative multi-agent game under uncertainty. Sufficient conditions for the existence of solutions of the two-stage SVI are provided. The numerical implementation is constructed by proposing a regularized SAA method. Furthermore, we prove the convergence of the regularized SAA method as the regularization parameter tends to zero and the sample size tends to infinity. A real application of crude oil market equilibrium is studied using our approach. More specifically, the oligopolistic market is characterized by a two-stage stochastic production and supply planning game and the market share is achieved at the equilibrium. Numerical results are presented based on historical data over the last fifty years and the effectiveness of our two-stage SVI approach of market share is demonstrated.

Appendix A Appendix

Proof of Theorem 3.7.

Proof A.1**.**

By direct computation, we have that

[TABLE]

Then, we can rewrite problem (11) as below:

[TABLE]

for $j\in\mathcal{J}$. From the first complementarity condition in (37), we have $(y^{\epsilon})_{j}$ as follows:

[TABLE]

for $j\in\mathcal{J}$. Similarly, we can derive that

[TABLE]

for $j\in\mathcal{J}$. Note that $(y^{\epsilon})_{j}=0$ implies $(y^{\epsilon})_{j}=0\leq x_{j}$, and we have $(\lambda^{\epsilon})_{j}=0$. Then, based on (38) and (39), we have for all three cases:

[TABLE]

where

[TABLE]

It follows that,

[TABLE]

which verifies (13). For the remaining of the proof, let $j\in\mathcal{I}_{2}$, we have

[TABLE]

and thus

[TABLE]

Analogously, we can derive from

[TABLE]

that

[TABLE]

Combining that of (40) and (41), we obtain

[TABLE]

Therefore, we have

[TABLE]

Substituting (42) into (41), we have derivation

[TABLE]

Then, we get

[TABLE]

that is,

[TABLE]

This completes the proof. ∎

Proof of Proposition 3.11.

Proof A.2**.**

Let $\hat{z}=(\hat{y},\hat{\lambda})$ be any solution of LCP$(q,M)$ and we have the derivation:

[TABLE]

where the second inequality follows from the positive semidefiniteness of $M$. Then, we have

[TABLE]

which implies the boundedness of $\lambda^{\epsilon}$,

[TABLE]

It follows from (43) that any accumulation point of $\{\lambda^{\epsilon}\}$ as $\epsilon\downarrow 0$ is the least $l_{2}$-norm solution. Since $M$ is positive semidefinite, we know from [32, Theorem 5.6.2] that there is a unique least $l_{2}$-norm solution. On the other hand, we know from Proposition 3.1, for any fixed $(x,\xi)$, $\hat{y}$ is unique. Therefore, the limit of $z^{\epsilon}$ exists as $\epsilon\downarrow 0$ and converges to the least $l_{2}$-norm solution of LCP$(q,M)$.

Due to the existence of limit for $z^{\epsilon}$ as $\epsilon\downarrow 0$, (15) can be derived directly from (13). In what follows, we focus on deriving the expression (16). To this end, for each $j\in\mathcal{J}$, three cases are discussed:

[TABLE]

Case 1: If there exists a sequence $\{\epsilon_{k}\}_{k=1}^{\infty}$ converging to [math] such that (44) holds, we have

[TABLE]

Thus, $\left|(\lambda^{\epsilon_{k}})_{j}-\bar{\lambda}_{j}\right|=0$.

Case 2: If there exists a sequence $\{\epsilon_{k}\}_{k=1}^{\infty}$ converging to [math] such that (45) holds, we have an estimation

[TABLE]

Thus, $\left|(\lambda^{\epsilon_{k}})_{j}-\bar{\lambda}_{j}\right|=0$.

Case 3: If there exists a sequence $\{\epsilon_{k}\}_{k=1}^{\infty}$ converging to [math] such that (46) holds, we have

[TABLE]

Thus, we have

[TABLE]

Collectively, we know from Case 1, Case 2 and Case 3 that

[TABLE]

Furthermore, we have that $T^{\epsilon_{k}}-\bar{T}\geq 0$ always holds. For the purpose of arriving at a contradiction, we assume $T^{\epsilon_{k}}-\bar{T}<0$. Then (48) implies that

[TABLE]

Moreover, (47) and (49) induce

[TABLE]

respectively. Clearly, we have $T^{\epsilon_{k}}-\bar{T}\geq 0$, which contradicts our assumption. In addition, we have

[TABLE]

Then, it follows that

[TABLE]

where the third inequality follows Lemma 3.9 and the continuity of $T^{\epsilon}$ that

[TABLE]

To summarize, for each $j\in\mathcal{J}$, we always have

[TABLE]

Then, according to the definition of $l_{2}$-norm, for any given $x\in\mathbb{R}_{+}^{J}$ one can compute

[TABLE]

∎