Roberts and Reagans (2007) show that the more a particular wine is exposed to critical ratings, the steeper the relationship between ratings and prices is; which we interpret as a stronger relationship between prices and quality: following this logic, our corpus of study consists of top Bordeaux wines. The sources used for selecting these wines were the 1855 ranking, the Graves ranking, the Saint-Emilion ranking, and prices for Pomerol wines. As the winemaking process differs between red and white wines, considering that more of the former are included in the different rankings across Bordeaux, the corpus was restricted to red wines. Finally, to reduce the price variability due to very localized events such as hail, only vineyards above areas of 5 ha were included in the corpus. The resulting corpus includes 59 different Bordeaux red wines, which belong to four Appellations d’Origine Contrôlée, henceforth named appellations. The complete list can be found in the Supplementary Appendix Table S1, and it is represented in Figure 1.

The long-term auction prices for 59 Bordeaux red wines over 53 different vintages were collected from the auction website IDealwine ¹ , which provides an index of average prices from the last auction sales. The data collection starts from vintage 1960, at which price records are available for a large majority of wines, and it ends at vintage 2013. The wines with more than five missing entries since 1960 were removed from the lists, leaving a total of 39 wines (see Supplementary Appendix Table S1). These prices will be our representation of the bottle’s long-term prices since each price point is an actualized price for a vintage over 9 years old. Figure 2 displays the evolution of average bottle prices over all vineyards of the four appellations.

Following Jones and Davis (2000), each vintage was divided according to the different phenological events, which are the milestones of the grape’s development. Indeed, the same weather condition can have a different impact on the grapevine depending on its growth phase, which can be captured only when partitioning the calendar in a physiologically relevant way. The most important phenological events, given here with their code on the BBCH scale (Lancashire et al., 1991) are budbreak (BBCH 07), flowering (BBCH 65), véraison (BBCH 85, the onset of ripening, marked by the changing of color of the grapes). These events can occur earlier or later across different Bordeaux appellations or for different grapevine varieties (e.g. Merlot is generally earlier than Cabernet), so the records obtained mention approximate dates. Even though the harvest is not a phenological event, harvest dates have also been included in the calendar, because they mark the end of the climate’s impact on the grapes. These historical dates were compiled from the records of Château Latour and the University of Bordeaux, to establish an approximate calendar of phenological events spanning the timeframe 1960–2017 (Table 1 and Figure 3).

The phenological events of each vintage provide a natural partition of the growing season in specific intervals:

• Budbreak—flowering: from budbreak to the first flowers

These intervals will be used to aggregate the weather data collected in the next subsection.

The weather data was gathered from the SAFRAN reanalysis of Météo France (Vidal et al., 2010), available in 8 km grid points, with daily granularity since 1958. The vineyards of each appellation were assigned weather data from one grid point: Saint-Émilion (Lon: −0.14, Lat: 44.91), for the Saint-Emilion and Pomerol appellations, Pauillac (Lon: −0.767, Lat: 45.182) for Médoc, and Léognan (Lon: 0.542/Lat: 44.756) for Pessac-Léognan. The selected weather variables shown are of classical use in the literature as predictors of wine quality. Based on van Leeuwen and Darriet (2016) who showed a highly significant correlation between wine quality and water deficit, the variable Water Deficit (WD), is calculated based on the simplified grapevine transpiration formula (Riou et al., 1994): W D = k * E T 0 − P (1) with E T 0 the Penman evapotranspiration, P the precipitations, and k = 0.3 from budbreak to flowering and 0.6 after flowering (van Leeuwen et al., 2004).

These weather parameters are then averaged over the five phenological intervals defined in the previous subsection to yield. Five different predictor features each. For instance, P (Precipitations) provides five different features: P: budbreak—flowering (average daily precipitations from budbreak to the beginning of flowering), P: flowering (average daily precipitations during flowering), P: flowering—véraison, etc. The same goes with all parameters of Table 2, except for Frost days: this parameter is not averaged over five phenological intervals but summed only over the phenological interval Budbreak—flowering.

As an accuracy benchmark against which to compare our predictive model, we collect critical primeur grades. The most influential critic in recent decades was without a doubt Robert Parker. His primeur grades were collected for the vintages 1994 through 2013, for 19 vineyards (list in Supplementary Appendix Table S2). Hence, our models will be evaluated on vintages 1994 to 2013 (20 vintages). Provided the tasted wines were young and still poised to evolve, each grade was only given as an interval and updated later to a single grade: as the goal here is to provide a benchmark of predictive performance, a unique grade is created by averaging the lower and upper bounds.

On our corpus of wines, average Spearman correlations between ratings and 2022 prices exhibit an increase across vintages (Figure 4), which is confirmed by a Kendall-Tau test (trend: 0.179, p-value: 0.0002): this means that the correlation of prices with ratings is stronger for recent vintages. It could be that there is a recent trend of relying more on expert opinion. But this more likely suggests that ageing vintages progressively decorrelate in value from early ratings, which we attribute to the fact that wine reveals its true quality with age, as mentioned by Storchmann (2012). This impact on early prices, and ulterior partial decorrelation, tends to confirm the short-term self-fulfilling prophecy effect of ratings already described in the literature (Ali et al., 2008). This side effect arguably boosts the performance of critical grades for price prediction, making it harder for predictive models to compete.

Literature on hedonic wine price functions generally recommends using the logarithm specification for the price (Oczkowski, 1994; Schamel and Anderson, 2003), which is supported by the better correlation of this specification with critical ratings (Oczkowski and Doucouliagos, 2015), hinting that it has a linearly more consistent variation with quality: as a result, the log specification is used hereafter. As our goal is to build a predictive model of wine prices, the real evaluation setting must be reproduced, in which the model must predict a previously unseen data point (here, the logarithm of a vintage’s price), and can only access data points from previous years. This is a case of ex-ante model testing, where the testing timeframe is posterior to the training timeframe, which prevents the model from accessing future trends during its training (Aristodemou and Tietze, 2018). Previous literature has only used in-sample testing instead to present their results (Jones and Davis, 2000; Jones and Storchmann, 2001; Jones et al., 2005; Ashenfelter, 2008; Baciocco et al., 2014; Corsi and Ashenfelter, 2019), meaning that the model was tested on the same data that it was trained on: this explains important performance gaps between their models’ results in this paper and previous evaluations. In the rest of this study, for the prediction of vintage v, each model is trained on a timeframe starting with the first available vintage (1960) and ending at vintage v-1, as displayed in Figure 5A. Each model is tested on the 1994–2013 time range (20 vintages), thus implying different training windows: for instance, vintage 1994 is predicted by models trained with the data from the years 1960–1993, and vintage 2004 is predicted based on the years 1960–2003. Evaluation metrics used are the Mean Absolute Error (MAE), the coefficient of determination, noted R 2 , and the Spearman rank correlation noted ρ_S (definitions can be found in the Supplementary Appendix Table S3). It is important to note that the coefficient of determination can be negative for all models in the setting of out-sample prediction, as mentioned e.g. by Chicco et al. (2021).

This section aims at selecting, out of all weather variables built in Section 2.1, a set of the best predictors of wine quality. To this end, the available weather variables for all vineyards are normalized, then concatenated into one cross-vineyard ensemble. The same is done for prices, which allows to compute statistics between any weather variable and overall prices. The Pearson R correlation with prices is computed for all weather variables, and the strongest correlations in absolute value are displayed in Table 3.

WD: flowering - harvest has the clearest correlation with prices. This is coherent with previous research where moderate water deficit is often associated with high-quality wine (Fraga et al., 2013; van Leeuwen and Darriet, 2016; Alem et al., 2019; Fayolle et al., 2019) because it leads to an increased grape tannin and anthocyanin content in various varieties (Duteau et al., 1981; Matthews and Anderson, 1988; van Leeuwen et al., 2009; Blank et al., 2019), as well as increased sugar concentration in fruit (Castellarin et al., 2007; Zsófi et al., 2011). Thus we will include this weather variable in our predictor set. The next most strongly correlated feature is P: véraison - harvest, which we will not consider for inclusion in our predictor set because it is too strongly correlated with WD: flowering—harvest (Pearson correlation: −0.75).

DTR: véraison - harvest has a strong positive correlation with prices. A high diurnal temperature range has already been linked to high quality (Gladstones 1992; cited in Jones et al., 2005), because it would be a sign of both a high diurnal temperature (crucial for berry ripening), and cool night temperatures enabling the production of the secondary metabolites associated with high-quality flavors (Tonietto and Carbonneau, 2004). However, this assertion has not been supported yet by any data, with some studies even concluding the opposite: it is therefore still debated to our knowledge (de Rességuier et al., 2020), and would necessitate further investigation to be a proven point. Nonetheless, as in the present case, the variable seems to have a strong positive impact on prices, it could have another relationship with prices, e.g. high DTR could only imply that the nights are cool, which could have a positive influence on quality. It will thus be used as a predictor of prices.

Finally, the negative correlation of P: flowering with prices makes sense. Flowering precipitations have been documented to cause climatic coulure and thus reduce yield (Blank et al., 2019). Prior literature displays little evidence of a reduction of quality, probably because an assessment of this impact is made more difficult by the lack of accessible phenology calendar data, but from our discussions with vintners in the Bordeaux area, the occurrence of strong rain during the flowering period is a very bad signal for the vintage’s quality.

The three variables: WD: flowering - harvest, DTR: véraison - harvest, and P: flowering have low multicollinearity, and they match the literature: they will thus be used as a set of inputs for predictive modeling, named S1. Growing season temperature has not been retained as a predictor, contrary to Ashenfelter (2008), because in line with the findings of Almaraz (2015), it did not exhibit a strong correlation with prices in recent years.

We also consider a classical set of variables used in literature (Ashenfelter, 2008), which we name S2. The third set S3 includes the square of the average growing season temperature, to allow for a second-order impact of temperature and capture a potential bell-shaped answer of quality to temperature following Jones et al. (2005). Table 4 summarizes the sets of weather variables used for predictive modeling.

Owing to the very nature of an agricultural yield prediction problem, where input-output couples can only be obtained once per harvest, the data at hand is sparse. This increases the risk of overfitting, namely the phenomenon by which a flexible model would adapt too much to the dataset variance and be unable of generalization (Hawkins, 2004). The Ordinary Least Squares model (Goldberger, 1962), noted OLS, has no hyperparameter and reduced flexibility, which reduces the probability of overfitting the training data. But its main advantage resides in the clear explanation that it gives of the relationships between predictor variables and the output: this is probably why this model is used in the overwhelming majority of econometric models in the literature (Jones and Davis, 2000; Esteves and Manso Orgaz, 2001; Jones and Storchmann, 2001; Jones et al., 2005; Haeger and Storchmann, 2006; Ashenfelter, 2008; Corsi and Ashenfelter, 2019).

The environment of a grapevine cultivar is described through the french word Terroir as the intricate relationship between climate, soil, and production methods (Seguin, 1986). For Bordeaux Grands crus, although the soils and cultivars remain mostly unchanged, deep evolutions are ongoing in both climate and production methods.

The 2.5°C increase in average growing season temperature in Bordeaux between 1960 and 2017 (illustrated for Saint-Emilion in Supplementary Appendix Figure S1), undoubtedly has significant effects on the plants, for instance causing phenological events to happen earlier (van Leeuwen and Darriet, 2016). On the side of the vine-growing methods, literature lists a myriad of innovations (see Gutiérrez-Gamboa et al., 2021 for a recent review), which seem to also have had a tangible effect, for instance by allowing the vineyard to maintain high fruit and wine quality until now (Gambetta and Kurtural, 2021) despite previous predictions of decline (Jones et al., 2005; Hannah et al., 2013).

Due to these long-term changes in the grapevine’s growing conditions, Almaraz (2015) evidenced that the impact of average growing season temperature on wine quality in Bordeaux has evolved over the last decades. We extrapolate this finding to hypothesize that other weather variables also have an evolving effect.

Then according to this hypothesis, in order to model the impact of certain weather parameters on wine quality, models with time-invariant effect cannot perform on long vintage series. Thus we want to find a model that can adapt its coefficients for a time-varying impact. However, the scarcity of the data at hand also raises the necessity of always keeping a memory of the oldest data points. Local Least Squares (LLS) kernel regression solves this dilemma, while still presenting the desirable properties of explainability and adversity to overfitting presented above. The goal of this method is to improve linear regression by applying more weight in the training to temporally close data points. This intuition is represented in Figure 5.

The problem of wine quality modeling is expressed as an extension of Equation 1 for several weather parameters. The logarithm of price of vintage v , noted y v , is estimated by the multivariate local linear estimator: y v = α v + ∑ i ∈ [ 1 ,   m ] β i ,   v   * x i , v with { β i , t } i ∈ [ 1 ,   m ] the coefficients associated with each parameter.

In the setting of LLS regression, we calculate the coefficients α and β used in the prediction of vintage v as the solution to the following optimization problem, which is simply a least-squares loss function, with the key difference that each term is attributed a different weight. This weight is determined by a kernel function K that decreases with the distance v   –   t between the predicted vintage v and the training data point of vintage t , and this distance is rescaled by a bandwidth parameter h : α K ,   h ,   v ^ ,   { β K ,   h , v } ^ i ∈ [ 1 ,   m ] = a r g m i n α ∈ R n ,   β ∈ R n * m   ∑ t ( y t − α t − ∑ i ∈ [ 1 ,   m ] β i ,   t * x i , t ) 2 *   K ( v − t h )

To solve this equation, we write it in matricial form, with the intercept α K ,   h ,   v ^ and coefficients β K ,   h ,   v ^ integrated into a single matrix β K ,   h ,   v ^ : β K ,   h ,   v ^ = a r g m i n β ( Y v − β X v ) T W K ,   h ,   v ( Y v − β X v ) With Y v the column matrix containing all training prices, the matrix of weather parameters X = [ 1 x 1 ,   t min ⋯ x m ,   t min ⋮ ⋮ ⋱ ⋮ 1 x 1 ,   t max ⋯ x m ,   t max ] , the coefficient matrix β = [ α v min β 1 ,   t min ⋯ β m ,   t min ⋮ ⋮ ⋱ ⋮ α v max β 1 ,   t max ⋯ β m ,   t max ] , and W K ,   h ,   v the diagonal matrix of weights d i a g ( K ( v − t min h ) ,   K ( v − ( t min + 1 ) h ) , … ,   K ( v − t max h ) ) .

Ruppert et al. (1995) give the closed-form solution to this optimization problem: β K ,   h ,   v = ( X v T W K ,   h ,   v X v ) − 1 X v T W K ,   h ,   v Y v Which in turn yields the estimated price: f K ,   h ,   v = β K ,   h ,   v X v = ( X v T W K ,   h ,   v X v ) − 1 X v T W K ,   h ,   v Y v X v Note that in the above equation, W K ,   h ,   v depends on a scalar bandwidth parameter h , which must be chosen depending on the data available—and in our setting, the prediction of each vintage v has access to data up to vintage v − 1 . We follow the cross-validation criterion, introduced by Clark (1977), which outperforms competitors in the case of multivariate prediction on small datasets (Köhler et al., 2014) as is the case here. Thus, the bandwidth h K ,   v for prediction of vintage v with kernel function K is set as the minimiser of the leave-out square error: h K ,   v = a r g m i n h ∈ R   ∑ t   ( β K ,   h ,   t , − t X t − y t ) 2 With β K ,   h ,   t , − t the leave-out estimator for year t , calculated from all available years except year t to prevent overfitting.

Widely used machine learning models have also been implemented to represent a panel of commonly used methods as a baseline against which to compare the OLS and LLS predictions. A simple Decision Tree (Quinlan, 1986), Random Forest (Breiman, 2001), Gradient Boosting (Friedman, 2001), and Single Value Regressor (Vapnik, 1995). All these implementations use the Python Scikit-learn package (Pedregosa et al., 2011). These models have been selected because they can provide decent performance on small datasets (here, the worst case is the prediction of vintage 1994, where training uses the 34 vintages from 1960 to 1993). But none of these models have been used yet -to our knowledge-in an academic study on grape price or yield prediction, and we do not expect them to yield good performance.