## 1 Introduction  
In many optimization or variational inequality (VI) problems, some data can depend on parameters, and it is important to analyze the regularity of the solution with respect to those parameters. For a detailed account of some classical results on the topic of parametric optimization, the reader can refer to [1]. The investigation of the Lipschitz behavior of solutions to convex minimization problems, from the point of view of set valued analysis, has also been carried out in the influential paper [2], and recently extended to abstract quasi-equilibrium problems in [3]. On a more applied context, in the case of linear programming, sharp Lipschitz constants for feasible and optimal solutions of some perturbed linear programs were derived in [4, 5]. Another influential paper focusing on the estimate of Lipschitz constants for the solutions of a perturbed system of equalities and inequalities is [6], while in [7], the estimate of the Lipschitz constant is connected with of the role played by the norm of the inverse of a non-singular matrix in bounding the perturbation of the solution of a system of equations in terms of a right-hand side perturbation. The topic of parametric sensitivity analysis for variational inequalities has been developed by many authors in recent decades, and the excellent monograph [8] contains many related results. Let us notice, however, that while the number of papers dealing with local parametric differentiability for variational inequalities is very large, the investigation of parametric variational inequalities under the mere Lipschitz assumptions on data have received less attention, in particular with respect to the problem of estimating the global Lipschitz constant of the solution. The local Lipschitz continuity of the solution of a variational inequality with parametric constraint set was investigated in [9], while more recently, the authors [10] derived an estimate of the global Lipschitz constant of the unique solution of a parametrized variational inequality, by making an extensive use of Lagrange multipliers techniques. Specifically, the Lagrange multiplier computations were used to derive the contribution of the parametric variation in the constraint set to the corresponding variation in the solution of the variational inequality. Interestingly enough, similar results were obtained in [11] with purely geometric methods. In [12], the authors improved the above-mentioned estimate but realized that, due to the complexity of the problem, their estimates were also far from being optimal. In this paper, we derive new estimates for the Lipschitz constant of the solution of a parametric variational inequality, which considerably improve the previous ones. We then combine our results with a variant of the algorithm proposed in [13; 14] to provide an approximation, on a given interval, of a univariate Lipschitz function whose analytic expression is unknown, but can be computed in arbitrarily chosen sample points (such a function is often called a black box-function). It turns out that the same algorithm also provides an approximation of the integral of the black-box function on the interval under consideration. Thus, we can use our estimate of the Lipschitz constant of the solution, to approximate its mean value on a given interval. Our results are then applied to investigate a class of Nash and generalized Nash equilibrium problems on networks. The specific class of Nash equilibrium problems known as Network Games (or games played on networks) was investigated for the first time in the paper by [PERSON] et al. [15]. Players are modeled as nodes of a graph, and the social or economic relationship between any two players is represented by a link connecting the two corresponding nodes. Two players connected with a link are called neighbors.