EconBase
← Back to paper

The Information Content of Taster's Valuation in Tea Auctions of India

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

110,660 characters · 30 sections · 24 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

The Information Content of Taster's Valuation in Tea Auctions of India

abstractTea Auctions across India occur as an ascending open auction, conducted online. Before the auction, a sample of the tea lot is sent to potential bidders and a group of tea tasters. The seller’s reserve price is a confidential function of the tea taster’s valuation, which also possibly acts as a signal to the bidders. In this paper, we work with the dataset from a single tea auction house, J Thomas, of tea dust category, on 49 weeks in the time span of 2018-2019, with the following objectives in mind: \begin{itemize} • Objective classification of the various categories of tea dust (25) into a more manageable, and robust classification of the tea dust, based on source and grades. • Predict which tea lots would be sold in the auction market, and a model for the final price conditioned on sale. • To study the distribution of price and ratio of the sold tea auction lots. • Make a detailed analysis of the information obtained from the tea taster's valuation and its impact on the final auction price. \end{itemize} The model used has shown various promising results on cross-validation. The importance of valuation is firmly established through analysis of causal relationship between the valuation and the actual price. The authors hope that this study of the properties and the detailed analysis of the role played by the various factors, would be significant in the decision making process for the players of the auction game, pave the way to remove the manual interference in an attempt to automate the auction procedure, and improve tea quality in markets.
center[center omitted — 158 chars of source]

Introduction

Tea (Camellia sinensis) is a manufactured drink that is consumed across the world. The tea crop has rather specific agro-climatic requirements that are only available in tropical and subtropical climates. Tea production, therefore, is geographically limited to a few areas around the world and is highly sensitive to changes in growing conditions. Majority of the tea producing countries are located in the continent of Asia with China, India, Kenya, Sri Lanka and Vietnam being the top producers (in that order), accounting for around 78% of the world tea production and 73% of exports. World tea production is estimated at over 5 million tonnes in 2015, valued around Rs 1 trillion. It increased by 4 percent to 5.2 million tonnes in 2014 (see FAO1 for further details). China remains the largest tea producing country with an output of 2.1 million tonnes in 2014, accounting for more than 40 percent of the world total, while production in India, the second largest producer, remained flat at 1.2 million tonnes in 2014, contributing around 30% of would production.

The tea industry is one of the oldest organized industries in India with a large network of tea producers, retailers, distributors, auctioneers, exporters and packers. Interestingly, India is also the world's largest consumer of black tea with the domestic market consuming around 1,000 million kg of tea during 2016. India’s annual production of tea is around 1,200 million kgs and the market size is estimated to be approximately Rs 20,000 Crore. The Tea Industry in India also derives its importance by being one of the major foreign exchange earners and for playing a vital role towards employment generation as the industry is highly labour intensive. India exports around 225 million kg of tea and is the fourth largest exporter in the world with Russia being its largest importer. The annual value of tea exports from India is around \$800 million. The other major importers of Indian tea are Iran, UK, Pakistan and UAE.

Tea is heterogeneous both over season and region, and even intra-region. Varieties of tea: 90% of the tea produced in India is of CTC (Cut, Tear and Curl) variety followed by Orthodox (9%) and Green Tea (1%). CTC tea is largely graded as Broken Leaf, Dust and Fannings. Each of these grades has a dozen of sub-grades based on the size of the grain etc. Several factors influence the demand for tea, including the price and income variables, demographics such as age, education, occupation, and cultural background. Apart from consumption, other main drivers of international tea prices are trends and changes in per capita consumption, market access, the potential effects of pests and diseases on production, and changing dynamics between retailers, wholesalers and multinationals (Source: FAO1).

Tea demand is very price sensitive. Price elasticities for black tea vary between -0.32 and -0.80, which means that a 10 percent increase in black tea retail prices will lead to a decline in demand for black tea between 3.2 percent and 8 percent, according to FAO. The average weekly volatility of CTC grade prices in Siliguri is 7%. Though the prices seem to be volatile, the tea prices show a specific pattern. The prices are at a peak as the new crop arrivals begin, with an increase in supplies, the prices witness a decline, lasting till the end of the season.

Due to the variation of quality, tea within the same grades are sold at a very wide range of prices in auctions. For example on any given auction day, say Broken Pekoe (BP) variety would sell at a range starting from as low as Rs 60 to Rs 250 per kg, depending on the producer mark (brand), quality and demand, making standardization difficult. Moreover, a given quality of tea is not available throughout the year. \footnote{The details regarding the tea industry has been collected from FAO1, FAO2 and website FAO3.}

{\bf Overview of e-auctions of tea:} An e-auction is a primary marketing channel for selling tea to the highest bidder. The auction system serves two basic purposes. The first purpose is to facilitate price discovery by bringing the buyers and sellers to a common platform with broker’s intermediation. Buyers bid for lots of tea and each lot is sold to the winning bidder. The second purpose is that the auction system provides a guaranteed transaction protocol for the transaction. The transaction includes activities such as delivery of tea to the warehouse, storing, sampling, bidding and payment (see FAO4 for a discussion).

Since September 2016, the auctions are pan India. It means that a member registered with any tea trade association anywhere in India can directly participate in any e-auctions conducted by Tea Board. Earlier, the buyer registered with local tea association could only buy from the e-auctions taking place in respective centers. The tea auction system brings the buyers and sellers together, to determine the price through interactive competitive bidding on the basis of prior assessment of quality of tea. Manufactured tea is dispatched from various gardens/ estates to the auction centres for sale through the appointed auctioneers, on receipt of which, the warehouse keeper sends an arrival a ‘weighment report’ showing the date of arrival and other details pertaining to the tea including any damage or short receipt from the carriers. The tea is catalogued on the basis of their arrival dates within the framework of the respective Tea Trade Associations, the quantities are determined according to the rate of arrivals at a particular auction centres. Registered buyers, representing both the domestic trade and exporters receive samples of each lot of teas catalogued, which is generally distributed a week ahead of each sale enabling the buyers to taste, inform their principals and receive their orders well in time for sale. The auctioneers taste and value the tea for sale and these valuations are released to the traders. Guidelines for the price levels likely to be established when the tea is sold are formulated on the basis of these valuations and last sale price. J. Thomas & Co. Pvt. Ltd is the largest auctioneer in the world, handling over 200 million kg of tea a year, which is one-third of all tea auctioned in India.

Given the above complexities, the report aims to evaluate the feasibility of automating the pricing process to the extent of dispensing with the manual testing and valuation steps. The next section describes the data set we have used for our analysis. We first discuss the clustering exercise according to Grade and source in Section 3. Section 4 discusses the pattern of volumes over different months of the year. The pattern of salability of different tea lots is investigated in section 5. Section 6 discusses the Price to value ratio to gain some insight into the pricing pattern which finally is used fully in the pricing models developed in Section 7. A detailed analysis of the value - price causality question is made in Section 8. The final comments on the feasibility of automation and future plans are discussed in section 9.

Description of data

In this report, we have used J-Thomas datasets on the weekly tea details for Kolkata Dust tea, Orthodox tea details, CTC tea details and Darjeeling tea details. Moreover, we have the e-auction statistics as a part and parcel of the dataset. We have used the data in 2018 for modelling, as the training data, and the 2019 data has been used for cross-validation. We begin our initial modelling assuming that the model, conditioned on the relevant factors, does not depend on the year of the auction. We shall later see in the cross-validation procedure, that our predictions are quite satisfactory to assert that our assumption was not falsified.

The e-auction statistics (2018-19) consists of the name of tea leaf type, total lots offered in auctions, total amount sold in packets and in quantity, and average price. For each such tea leaf type, detailed info on the weekly sale, total amount sold in packets and in quantity, and average price has also been provided. An initial overview of the characteristics of the data at hand, has been produced in our_earlier.

In the J-Thomas datasets, we find lot numbers (hence the difference between the maximum and minimum lot number would give us the number of lots offered), the categorical variable: the grade of the tea, number of packages offered, the valuation given by the agency, and finally the auction selling price.

In our dataset, we have 25 types of tea grades available names, namely,

center[center omitted — 954 chars of source]
comment[1] "CD" "CD1" "CHD" "CHD1" "CHU" [6] "D" "D(F)" "D(SPL)" "D1" "D1(SPL)" [11] "GTDUST" "OCD" "OD" "OD(S)" "OD1" [16] "OPD" "OPD(CLONAL)" "OPD1" "ORD" "PD" [21] "PD(FINE)" "PD(SPL)" "PD1" "PD1(SPL)" "RD1"

However, there are only three main broad categories for tea dust: D1, PD and PD1 according to Wikipedia wiki. Hence our grades require clustering. Out of the 49 weeks of data we possess, 38 weeks of auctions are from 2018, from the period of January to December. However, since not every week had an auction, we thus have data on weeks 2-9, 18, 20, 22-38, 40-41,43-52. Along with this, the remaining data is on the first 12 weeks in 2019, (except the 11th week, in which no auction took place), which we shall be using for cross-validation. Hence, based on our training set, we club together the tea grades for which there are less than 38 observations in the entire dataset (as we require at least one observation per week on the average.) This clubbing is done by its proximity to the other tea grades, where the proximity is based on the qualitative similarity of the grades. Thus we obtain the following 14 clubbed grades till now.

itemize• CD: Consisting of CD, CHD and CHU • CD1: CD1 and CHD1 • D: D and D special • D(FINE) • D1: D1 amd D1 Special • OCD • \textbf{OD}: OD and OD-special • \textbf{OD1} • \textbf{OPD}:OPD, OPD-Clonal and ORD • \textbf{OPD1} • \textbf{PD}: PD and PD-Special • \textbf{PD(FINE)} • \textbf{PD1}: PD1 and PD1-Special • \textbf{RD1}

Due to the only packet of GT Dust in the dataset, and that too remaining unsold, and further, due to its lack of immediate similarity from any of the existing tea grades, we conclude that GT Dust is a very rare category in this dataset, and hence we leave it out for the classification problem for now.

commentTemporally, we are presented with the data of the following weeks in 2018 (in the form of week numbers in which they occur in a year): $$2 \ 3 \ 4 \ 6 \ 7 \ 8 \ 9 \ 18 \ 20 \ 22 \ 23 \ 24 \ 25 \ 26 \ 27 \ 28 \ 29 \ 30 \ 31 \ 32 \ 33 \ 34 \ 35 \ 36 \ 37 \ 38 \ 40 \ 41 \ 43 \ 44 \ 45 \ 46 \ 47 \ 48 \ 49 \ 50 \ 51 \ 52$$ This has been the basis for the model, put into brackets of months, and has been discussed later. For 2019, we have the data on the following weeks, presented in the same format as that in 2018: $$ 1\ 2\ 3\ 4\ 5\ 6\ 7\ 8\ 9\ 10\ 12$$ These have been used as our cross-validation data. The analysis can be found in the subsequent sections.

Classification and Clustering

Clustering based on Grade

The number of clusters so formed in the previous subsection is still quite large, and we suspect that they have quite an inherent similarity in their characteristics and hence in their market appeal and corresponding auction transactions. To have an idea about this, we form a dissimilarity matrix among the grades to visualize the measure of degree of similarity across the clusters. We form the dissimilarity matrix based on the Volume Weighted Valuations and use the metric: $$\text{Dissimilarity} (d)=2(1-\rho^2)$$ where $\rho$ is the product moment correlation correlation coefficient between these time series of Volume Weighted median Valuations for different grades. Thus we create the dissimilarity matrices between these grades, and obtain Figure (ref), where the darker shade represents a larger value of dissimilarity, while a lighter shade shows that the clusters are similar.

figure[figure omitted — 252 chars of source]

One approach of clustering these is using hierarchical clustering based on this dissimilarity measure, to obtain the following clustering dendrogram in Figure (ref).

figure[figure omitted — 217 chars of source]

Analogously, considering Volume weighted mean valuations, we obtain Figure (ref).

figure[figure omitted — 222 chars of source]

However, note that correlation is invariant to change of scale and origin, hence if a particular grade has even twice a valuation than another, using correlation as a measure of clustering would essentially nullify the effect of even twice the valuation, and would land the two grades into the same cluster, in spite of the fact that the grades are quite distinct in their market characteristics.

Furthermore, we need to incorporate the idea that we have the data for multiple weeks, and the clustering should involve the clubbing for all the weeks combined.

Thus we use the following idea:

itemize• First we use Bayesian Information Criteria to figure out the number of clusters so that the model has the largest information. This figure came out to be $6$. • Then, we used a Gaussian Mixture Model and ran the Expectation-Maximization Algorithm (often known in the literature as the EM-GMM clustering). We performed clustering both based on median as well as the mean, these two resulting in slightly different clusterings. We stick with the median based clustering due to its robustness. • Both the price and the valuation was considered when applying EM-GMM clustering, in order to effectively capture the whole pattern present in the data.

Then, we define a new similarity structure, where the $i$-th and the $j$-th grade's similarity is proportional to the number of weeks they have occurred in the same cluster by EM-GMM method. Thus, we form a similarity matrix with $(i,j)$-th entry of the matrix characterizing the similarity in terms of number of co-occurrence weeks. The Mosaic plot corresponding to these similarity matrices is given in Figure (ref) and Figure (ref). Finally based on this similarity matrix, we conduct a hierarchical clustering to obtain the clusters, shown in Figure (ref) and Figure (ref).

figure[figure omitted — 315 chars of source]
figure[figure omitted — 320 chars of source]
figure[figure omitted — 209 chars of source]
figure[figure omitted — 215 chars of source]

Thus we obtain the following 6 clusters based on grade:

itemize• Cluster 1: OD, OD-Special, OPD1 (3 grades) • Cluster 2: OCD, OCD1, OD1 (3 grades) • Cluster 3: D-Fine, CD1, CHD1, RD1 (4 grades) • Cluster 4: D, D-Special, CD, CHD, CHU, PD, PD-Special (7 grades) • Cluster 5: PD-Fine (1 grade) • Cluster 6: OPD, OPD-Clonal, ORD, D1, D1-Special, PD1, PD1-Special (7 grades)

while the GT Dust category has been left out due to lack of sufficient data points.

Diagnostics of Grade Clustering

As mentioned before, according to Wikipedia wiki, the three main categories for dust type tea leaf are D1, PD and PD1. However, our clustering algorithm shows that Tea Grade D1 and PD1 are similar in market characteristics by putting them into same cluster. On this note, we consider two different clusterings which might be possible.

enumerate• Method 1: The 6 clusters as mentioned above. • Method 2: 8 clusters, cluster 1, 2, 3, 5 remaining as it is, while cluster 3 gets broken into two separate clusters, one containing grades PD, PD-Special and another containing CD, CHD, CHU, D, D-Special. Similarly, we divide cluster 6 into two separate clusters, one containing OPD, OPD-Clonal, ORD, D1, D1-Special and another containing PD1, PD1-Special.

To find out which one of the above clusterings would be better for further analysis, we find out the proportion of total variation of both weekly price and weekly valuation which is explained by the clusters. Some of the obtained results for both of the clusterings are given in Table (ref). It was found that, for both weekly price and weekly valuation, about 50% variation of these variables over different lots are explained by the clusters alone. Also, using 8 clusters instead of 6 clusters increases this proportion of explained variation by at most 2%, which is not significant in contrast to the loss of simplicity of subsequent works. This diagnostic suggests us to stick with the original clustering, with 6 clusters as defined previously.

table[table omitted — 779 chars of source]

Clustering by Source

Now, note that, if we just cluster the tea based on their grades, then we are foregoing valuable information about the source from which the tea has been produced. This information would be very relevant for our subsequent discussions, and it would not be a good idea to get rid of it. Valuations about a tea grade depend on gardens they are originating from, as it utilizes the idea about the environment used for their nourishment, soil levels, and other significant factors. Hence the source of the tea dust is of utmost importance. However, there are 238 tea gardens (293 including their Clonal, Royal, Gold and Special variants) from which the tea has originated, and again as above, we suspect that maintaining track of each of the tea gardens would be intractable as well as redundant, as the tea also show similarity in characteristics, a significant factor of which might be geographical proximity. Hence we undergo clustering based on the volume-weighted median to obtain dendrogram given in Figure (ref). These dendrogram has again been clustered using the EM-GMM algorithm and then the time-based similarity matrix as in the preceding section.

figure[figure omitted — 201 chars of source]
figure[figure omitted — 264 chars of source]

Before going into the final source-based clusterings, we need to keep the following things in mind as well:

itemize• The districts for the source of the tea dust vary from the northern fringes of West Bengal to the entirety of Assam. Hence it might be a good idea to keep the districts of West Bengal and Assam well segregated. • Geographical proximity might be a concern for the explanation of the districts that belong to the same cluster. • Topography, soil structure, and water source may also be the reason for similar or dissimilar market characteristics of the tea dust. • Most of the tea dust from West Bengal comes from Jalpaiguri, as portrayed in the data, while the other districts have a significantly lesser number of such packets to be sold. Hence it would be good if the districts of West Bengal are classified keeping this in mind.

We attach here the maps of the tea-dust producing districts of West Bengal and Assam to give a visualization of the geographical proximities (in Figure (ref) and Figure (ref) ).

figure[figure omitted — 201 chars of source]
figure[figure omitted — 200 chars of source]

Thus we finally use the following 7 clusters based on the source of the tea dust.

itemize• Cluster 1: Darjeeling, Cooch Behar (Koch Behar), Uttar Dinajpur, Jalpaiguri (4 West Bengal districts) • Cluster 2: Karimganj, Hailakandi • Cluster 3: Bongaigaon, Cachar, Udalguri, Darrang, Dima Hasao • Cluster 4: Lakhimpur • Cluster 5: Nagaon • Cluster 6: Sivasagar, Tinsukia, Golaghat, Jorhat • \textbf{Cluster 7}: Baksa, Dibrugarh, Sonitpur

Distribution of Volume over Different Months

Once the clustering for tea grades and sources are obtained, the next idea would be to create a temporal clustering. The data for 2018, presented in the form of weeks, are put into buckets of months. This is done to capture the seasonal variations in the market characteristics of tea dust. But again the exact week number might be too redundant an information, as the tea dust appearing in the market does not fluctuate as frequently as weeks, but might vary by seasons. Thus, to incorporate such a possibility of market dynamics, we calibrate the data of tea dust grades by months. This gives us the Table (ref). To obtain Table (ref), we find out the number of tea packets offered in a lot and multiply it with the net average weight of the tea packets to obtain the total amount (or volume) of tea offered. Then, for each cluster of grade, we find the proportion of its total volume which is offered during a specified month. This gives us a basis for clustering to analyze the supply side of market dynamics.

table[table omitted — 1,279 chars of source]

The mosaic plot for the above table has been included in Figure (ref) for a better visualization.

figure[figure omitted — 267 chars of source]

Predicting Salability of Offered Tea lots

For analyzing the auction of tea packets, we should concern ourselves with the proportion of tea packets to be sold, and relate its valuation and several other characteristics to it. A successful attempt at predicting the probability of being sold (or being unsold) of an incoming tea packet based on its Grade, Source, Valuation and the current month, would give an insight for automating the auction process, alongside enabling an opportunity to study the effect of Valuation on determining the market characteristics of tea grades.

A Primary Inspection

Primary inspection is made to see whether any particular type of tea grades are more likely to be sold at the auction than other types. Table (ref) shows the proportion of tea lots being sold and the total number of tea lots offered across different grades.

table[table omitted — 1,069 chars of source]

Note that, Table (ref) shows that Fine variant of a tea grade is offered more rarely than its original variant, and its probability of getting sold at the auction also increases. Also, the Special variant of any tea grade is rarer to be offered than its Fine variant, and for this reason, there is not a significant number of observations to conclude whether it increases or decreases the probability of getting sold at the auction.

To find out whether the probability of being sold significantly depends on the variant of tea grades, we perform a simple one-way Analysis of Variance model with the indicator of being sold as the response variable and the type of variant as our treatment variable. We obtain the results as shown in Table (ref). Clearly, the probability of being sold for Fine variants of tea grades are higher than regular variant, and the p-value is small indicating that there are a significant number of observations to support this. However, the proportion of being sold is possibly lower in Special variants than in Regular ones, but higher p-value indicates that there is not enough evidence to support this.

table[table omitted — 495 chars of source]

Similar to this, we tried to find out whether different variants of Source (for example, Clonal, Gold, Royal, etc.) affect the probability of getting sold. Again we perform an Analysis of Variance model, however, with the variant of the tea garden (or Source) as our treatment effect. We obtain the results as shown in Table (ref). From this, we note that if the tea packet has come from a Clonal tea garden, its selling probability is expected to be higher than Regular ones by 0.044, and the smaller value of p-value indicates evidence to support this claim. Similarly, the tea packets produced from the Gold type variant of Garden is expected to be 15% less probable to be sold at the auction. Also, with 95% confidence, we can say that the tea packets produced from the Royal type variant of Garden are 22% less likely to be sold.

table[table omitted — 686 chars of source]

Model for Prediction

To build a predictive model to predict whether a tea packet will be sold at the auction or not, based on its Valuation, Grade, Source, and the current time, we use three competing models.

enumerate• Logistic Regression • Generalized Additive Model • Mixture of Logistic Regression

We divide the total dataset of the year 2018 into training and cross-validation sets, with the training set containing 70% of the samples. The cross-validation set is used to select the model. The dataset of the year 2019 is kept as a testing set, which is used to evaluate the performance of the finally selected model for prediction. The percentage of sold tea lots is kept almost similar for both the training and testing sets about 81%. For each of the sets, we consider each combination of Valuation, Grade, Source, and Month, and obtain the proportion of sold tea lots among all tea lots offered under that combination. On the other hand, the predicted probability of being sold under that combination is estimated from the trained model. Both of these probabilities are visually and analytically compared against each other to assess the performance of the trained model for both sets. For analytical comparison, we use three measures as follows;

enumerate• Null and Residual Deviance: For a generalized linear model, $$\underbrace{\text{Null deviance}}_{\text{df}}=\underbrace{2\left(\Lambda(\text{Saturated Model})-\Lambda(\text{Null Model})\right)}_{\text{df(Saturated Model) - df(Null Model) }}$$ and $$\underbrace{\text{Residual deviance}}_{\text{df}}=\underbrace{2\left(\Lambda(\text{Saturated Model})-\Lambda(\text{Proposed Model})\right)}_{\text{df(Saturated Model) - df(Proposed Model) }}$$ where $\Lambda(\cdot)$ stands for the log-likelihood under the model in the argument. The saturated model is characterized by each data point bringing in its new parameter, ie, it provides no further scope for the addition of parameters corresponding to the data points. Thus all $n$ parameters are to be estimated. The null model is the exact antipodal in the sense that it assumes a common parameter for all the data points, and thus only one parameter needs to be estimated. The proposed model is the model one proposes, with $p+1$ parameters. A small null or residual deviance indicates that the corresponding one or p+1 parameter model explains the data well. Formally, under a truly good fit, the deviances follow a chi-squared distribution with the mentioned degrees of freedom. • RMSE: Root Mean Squared Error measures the $L_2$ distance between the predicted probabilities and the true probabilities. Hence smaller values are preferred. • MAE: Mean Absolute Error measures the $L_1$ distance between the predicted and true probabilities. Hence smaller values are preferred.

Logistic Regression

We began with attempting to fit a logistic regression model with the full data. The following results came up:

itemize• Null deviance: 18055 on 18381 degrees of freedom. • Residual deviance: 17813 on 18359 degrees of freedom. • Null deviance - residual deviance = 242 with degrees of freedom 22. • Pseudo $R^2$: 0.01340349. • RMSE is 0.295 for training set and 0.312 for cross-validation set. • MAE is 0.214 for training set and 0.235 for cross-validation set. • Figure (ref) shows how bad the fit is for training and testing sets respectively.
figure[figure omitted — 220 chars of source]

Hence this model fails miserably in predicting whether a packet would be sold.

Generalized Additive Model

We use a Generalized Additive Model with a binomial family, with a smooth cubic spline fitted on the Valuation of tea packets as a predictor. The following results came up:

itemize• Null deviance: 18055 on 18381 degrees of freedom. • Residual deviance: 17627 on 18357 degrees of freedom. • Null deviance - residual deviance = 429 with degrees of freedom 24. • Pseudo $R^2$: 0.02373118. • RMSE is 0.287 for the training set and 0.306 for the cross-validation set. • MAE is 0.21 for the training set and 0.2322 for the cross-validation set. • Figure (ref) shows some improvement over logistic regression, but still, the fit remains too bad to be of any use.
figure[figure omitted — 220 chars of source]

Mixture of Logistic Regression

We try using a mixture of logistic regressions to predict the selling potential of the packets. The model, in general, for a mixture of $S$ components, is given by logit_mix (and extended in logit_mix_2, logit_mix_3) \[ H(y|T,\mathbf{x},\mathbf{w},\mathbf{\Theta})=\sum_{s=1}^S \pi_s(\mathbf{w},\mathbf{\alpha})\text{Bi}(y|T,\theta_s(x)) \] where $\mathbf{w}$ stands for the concomitant variables on which the mixing proportions $\pi_s$ depend, $\text{Bi}(y|T,\theta_s(\textbf{x}))$ is the binomial distribution with number of trials equal to $T$ and success probability $\theta_s\in(0,1)$, is modelled by usual logistic modelling; $\text{logit} \left(\theta_s(\mathbf{x})\right)=\mathbf{x}^\top \mathbf{\beta}^s$. The concomitant variable is assumed to have a multinomial logit model, i.e. of the form $$\pi_s(\mathbf{w},\mathbf{\alpha})=\dfrac{e^{\mathbf{w}^\top\mathbf{\alpha_s}}}{\sum_{u=1}^S e^{\mathbf{w}^\top\mathbf{\alpha_u}}} \ \forall s$$

Here, we have used our independent variable $x$ to be the corresponding valuations, $T$ being the number of packets arriving, and success denoting the event that the packet is sold. The concomitant variables are the source, grade cluster and month of the packets.

Here we consider 2 to 5 component mixture for this. Based on the Bayesian Information Criterion, the 3 component mixture of logistic regression is chosen which yields a BIC value 6183.014.

itemize• RMSE is 0.222 for the training set and 0.2488 for the cross-validation set. • MAE is 0.1522 for the training set and 0.1811 for the cross-validation set. • Figure (ref) provides the plots for fits, which shows significant improvement over the previous models. Figure (ref) gives the rootogram of the components, and the peaks at 0 and 1 suggest how well these components can be identified. Note that, the observations assigned to component 1 are marked in the rootogram. The peak at 1 for component 1 shows that this component is well separated from others. • Taking into consideration the metrics for fit as well as the plots, we finalize this model to predict the probability of a particular packet of tea being sold, given its corresponding covariates. • As we decided to use this as our final prediction model, we can update its parameters based on all of the training and cross-validation set. Finally, we evaluate its performance on the testing set (comprises of data in 2019). \begin{itemize} • Updated model achieves RMSE of 0.214 and 0.2519 in the whole training set and testing set respectively. • This achieves an MAE of 0.146 and 0.1684 in the whole training set and testing set respectively. • Figure (ref) provides the plots for the fitted model for both training and testing sets. • Table (ref) provides the summary of this model. \end{itemize}
figure[figure omitted — 265 chars of source]
figure[figure omitted — 214 chars of source]
figure[figure omitted — 287 chars of source]
table[table omitted — 961 chars of source]

Distribution of price to valuation ratio for grade clusters

Valuation by experts does provide a significant knowledge of what the final price of the transaction would be. The entire process runs with the base price being set at some fixed proportion of the valuations, and thus the following price of transaction revolves significantly across this measure. In this section, we would attempt to fit distributions over the ratio of price and valuations to account for its variability and shape of distribution curves. We have attempted this exercise with the natural logarithm of the ratio of price and volume, to have full support over the real numbers. Observe that if, $\log\left(\dfrac{X}{Y}\right)\sim N(\mu,\sigma^2)$, then $$\mathbb{E}\left(\dfrac{X}{Y}\right)=\exp\left(\mu+\dfrac{\sigma^2}{2}\right)$$ and $$\mathbb{V}\text{ar}\left(\dfrac{X}{Y}\right)=\exp\left(2\mu+\sigma^2\right)\left(e^{\sigma^2}-1\right)$$ The histograms of the ratio of price and valuations for various grade clusters as obtained from the data are given in Figure (ref) and Figure (ref).

figure[figure omitted — 211 chars of source]
figure[figure omitted — 212 chars of source]

For Cluster 1 (OD, OD-Special, OPD1)

We have attempted fitting a single normal distribution as suggested by the histogram. However, this results in a poor fitting of the data. Hence, we attempt to fit a mixture of two normal distributions to the data. The chi-squared goodness of fit yielded a p-value of 0.2274 and the Kolmogorov-Smirnov statistic yielded a p-value of 0.1588, which is quite satisfactory for our case. Thus we report the distribution of the ratio of price and valuation of Cluster 1 to be a mixture of two log-normal distributions with the following properties

center[center omitted — 281 chars of source]

For Cluster 2 (OCD, OCD1, OD1)

The histogram shows a single modal distribution. We tried fitting a single lognormal distribution, our p-values for the chi-square goodness of fit statistic came out as 0.2665 and for One sample Kolmogorov Smirnov test came out as 0.6908, which suggests a reasonably good fit. Thus we report the ratio for this cluster to follow a an unimodal log-normal distribution with the following properties.

center[center omitted — 200 chars of source]

For Cluster 3 (D-Fine, CD1, CHD1, RD1)

The same pattern follows as in Cluster 2, unimodal from histogram, but a single log-normal fits badly (p-value 2e-04). But fit with mixture of two log-normals give reasonably well fits (p-value 0.3855). Thus we report a mixture of two log-normal distributions with the following properties.

center[center omitted — 295 chars of source]

For Cluster 4 (D, D-Special, CD, CHD etc.)

The histogram shows a unimodal and highly leptokurtic structure, thus we attempt to fit a single log-normal distribution. However, it fits badly. A mixture of two log-normal distributions fits reasonably better, yielding a p-value of 0.4998 for Pearsonian goodness of fit test, and a p-value of 0.324 for Kolmogorov Smirnov test. Hence we report a mixture of two log-normal distributions with the following properties. The closeness of the means and smaller variances in both the components explain the reason for a unimodal looking structure in the histogram.

center[center omitted — 294 chars of source]

For Cluster 5 (PD-Fine)

In this case, as the histogram shows a bimodal shape, we try fitting with a 2 component mixture of log-normal distribution. It gives reasonably well p-value, 0.8231 for Pearson's chi-sqaured goodness of fit test and 0.9891 for Kolmogorov Smirnov's test. Hence we report the distribution of the ratio to be a mixture of two log-normal distributions with the following parameters.

center[center omitted — 297 chars of source]

For Cluster 6 (OPD, OPD-Clonal, ORD, D1 etc.)

The histogram for this cluster is characterized by its unimodality and slight positive skewness. A single log-normal yielded a p-value of 0.0104. Hence, we moved on to a mixture of two log-normal distributions, and this time, the p-value came out to be a higher value of 0.06225. We also tried to fit a three-component mixture of log-normal distributions, however, that does not increase the p-value by a significant amount. Thus we report a mixture of two log-normal distributions with the following properties, to maintain the simplicity of the underlying model.

center[center omitted — 299 chars of source]

Analysis of the Price Model

To model the pricing system of the tea market, we consider modeling the demand side by consideration of the Valuation of tea grades, its grade, source, and the month in which the tea lots are available. On the other hand, to model the supply side of the market, we consider the volume of the tea lots as our main predictor. Therefore, our pricing model should include these variables.

To check whether the variant of the tea gardens (Clonal, Gold, Royal, Special, etc.) should be included in the pricing model, we simply fit a one-way Analysis of Variance model with Price as our response variable and the variant of tea garden as the possible treatment variable. It is found that this factor explains a sum of squares of 842405 with 4 degrees of freedom, yielding an F-statistic value of 102.41 and consequently extremely small p-value. Therefore, based on the data, we find sufficient evidence to incorporate this factor into our pricing model in order to have better predictability.

Nonparametric Exploration of limitations of a Pricing Model

Before specifying a statistical model for the prediction of Valuation and Price of a tea packet based on several of its characteristics, it is extremely important to understand the limitations of how much we can do first. To this end, it is found that there may be several tea packets with exactly the same characteristics, (coming from the same garden, is of the same grade and comes in the same week), and between them, the Range of their Valuations and Prices are calculated. If we consider the empirical CDF of such all possible ranges, then as seen from Figure (ref), in most cases, those packets are subjected to exactly the same Valuation by the auctioneers, however, the final Prices at they are sold may be very different. This exploration suggests that if we simply do away with Valuation and output a single prediction of Price for tea packets based on its Grade, Source, and Week of the year when it is held for the auction, we can almost hope for $55\%$ accuracy in prediction. However, if we wish to have a $90\%$ accuracy in predicting price, we must allow approximately $30$ rupees of deviation from the actual price. Even with a more robust measure of dispersion, mean deviation about median, the conclusion of this exploration remains the same, as seen from Figure (ref).

figure[figure omitted — 352 chars of source]
figure[figure omitted — 326 chars of source]

Since there are lots of Grades, Source and possible weeks combination provided in the dataset, modeling a different prediction for each of these combinations would require a great number of parameters in the model. However, as we have already performed a reasonable clustering analysis of different tea grades and the source gardens, we wish to explore whether having a prediction for tea grades with similar cluster characteristics will be within the economical tolerance level for the auctioneers. Figure (ref) and (ref) shows the corresponding ECDFs of Ranges and Mean deviations about median of Valuations and Prices, for tea packets sharing same cluster characteristics. As seen from them, such a model that predicts at the cluster level would be inadequate in modeling either the valuation or the price soundly.

figure[figure omitted — 384 chars of source]
figure[figure omitted — 358 chars of source]

Linear Price Model

For ease of interpretability, we start with a simple linear model for pricing system with the grade, source clusters, month of availability, variant of source garden, the volume of the tea packet, and valuation of the tea packets as our predictors. The model is given by;

comment\begin{multline*} Price=\beta_0+\underbrace{\beta_1Grade+\beta_2Source Clusters+\beta_3Month+\beta_4Garden Variant+\beta_6Valuation}_{\text{Factors for Demand}}\\+\underbrace{\beta_5\text{Volume}}_{\text{Supply}}+\varepsilon\\ \text{ where } \varepsilon\sim N(0,\sigma^2) \text{ independently and identically distributed} \end{multline*}
align*[align* omitted — 358 chars of source]

The results from this model are summarized in Table (ref).

table[table omitted — 781 chars of source]
figure[figure omitted — 291 chars of source]
enumerate• Multiple R-squared for the fitted model is 0.9234, suggesting a very strong linear relationship among the predictors. • The Analysis of Variance decomposition between different effects has been shown in Table (ref). Clearly, each of the predictor variables explains a lot of variation in the price, and the extremely small p-values indicate that each one of the variables' contribution is significant. • The true price and predicted price based on the fitted linear model has been shown in Figure (ref). Note that, the predicted prices are randomly dispersed on both sides of the reference line, thereby supporting the assumption of homoscedasticity. Also, the standardized residuals have quantiles similar to that of the theoretical quantiles of a standard normal distribution, as assumed by the model, other than some drastic outliers present in both ends. • Based on the evaluation of the pricing model in the testing dataset, the 2.5% quantile of the residuals is -18.59145, while the 97.5% quantile of the residuals is 23.43802. Hence, this pricing model approximately makes an error about 20 Rupees, to both positive and negative sides, considering a robust measure of variation. • On the other hand, considering a classical approach to measure the standard error, we find the interval, with predicted value as the center and an error of 24.3229 added (or subtracted) to both sides contains the true price for 95% of the time.

Comparatively, proceeding with our main objective, we remove the valuation as predictor to see how much it affects our original linear pricing model. In this case, we find that a simple linear model with Price of the tea packets as response variable would make the residuals to be heteroscedastic. Therefore, we apply a variance stabilizing logarithmic transformation and use the natural logarithm of price of tea packets as response variable. Thus the model is given by;

align*[align* omitted — 248 chars of source]

where $\varepsilon\sim N(0,\sigma^2)$ independently and identically distributed. The results obtained are summarized in Table (ref)

table[table omitted — 1,856 chars of source]
figure[figure omitted — 299 chars of source]
enumerate• Multiple R-squared for the fitted model is 0.6431, suggesting a moderately strong linear relationship among the predictors. • The Analysis of Variance decomposition between different effects has been shown in Table (ref). Clearly, each of the predictor variables explains a lot of variation in the price, and the extremely small p-values indicate that each one of the variables' contribution is significant. However, to our surprise, the effect of Volume has been explained by other variables to some extent, although the volume of tea packets generates the supply side of the market. • The true price and predicted price based on the fitted log-linear model has been shown in Figure (ref). Note that, the amount of error in prediction increases as the true price increases, thereby supporting the evidence of heteroscedasticity as previously mentioned. Also, from the Q-Q plot, it is evident that there are some possible outliers present at lower prices. • Testing the performance of the fitted model on the testing set, the 2.5% quantile of the residuals comes as -0.29539, while the 97.5% quantile of the residuals is 0.37754 in logarithmic scale. Hence, this pricing model approximately makes an error from 74.43% to 145.86% of the true prices of the tea packets. This is based on a robust approach to estimate the standard error. • On the other hand, considering a classical approach to measure the standard error, we find that the predicted prices lie between 71.09% and 140.65% of the true prices about 95% of the time.

Thus the valuation part is significant for the prediction of the prices, and the process cannot be automated by a simple linear or log-linear model approach.

Causal Analysis between Price and Valuation

From the simple statistical analysis with a linear model performed above, it seems Valuation is indeed a pertinent variable in the explanation of the Pricing system. However, we could not yet allege that the auctioneers' valuation had a causal impact on the price level, although it seemed to be an indispensable predictor. To evaluate the causal impact of valuation on price, we lay down a part of the causal graph structure in Figure (ref). Note that there may be some causation structure between the variables (for example, not all grades occur in all months, hence there should be an arrow from month to grade). Nevertheless, the variables of interest are the valuation and the final price, and the parents to them are known, assuming we subscribe to a linear understanding of time and causality. Thus all the other variables in the model are parents to both valuation and price. Furthermore, valuation could be a parent of price. But, if the valuation is indeed an 'educated guess' of price as was originally intended, then, since the guess is based on only these variables; conditioned on these variables, price and valuation should be independent trials from (possibly same) distribution. Hence, they should be independent, had valuation not been a cause of price. Hence we aim to test this independence.

Causal Analysis with linear effect model

figure[figure omitted — 193 chars of source]

Thus, we wish to test whether $$\text{Valuation}\perp\!\!\!\perp \text{Price} \mid (\underbrace{\text{ Source, Grade, Volume, Garden, Month}}_{\text{rest}}) ?$$ However, we are faced with the fact that conditioning on so many variables render fewer data to provide reliable estimates for any inference. Hence, we take a different strategy, which is generally often taken when the conditioning variables are continuous. We model the log of valuations by a linear function of the other variables, and similarly for the log of price. That is

align*[align* omitted — 323 chars of source]

If these models provide a good fit, then we can look at the correlation between the residuals to identify the presence of a causal link between valuation and price. This is because the residuals in both the models are independent of the rest, and for two variables $A, B$ and $C$, we have, $$A\perp\!\!\!\perp B |C \Leftrightarrow [A - f(C)]\perp\!\!\!\perp [B - g(C)] | C \Leftrightarrow A - f(C) \perp\!\!\!\perp B - g(C) $$ if $A-f(C)$ and $B-g(C)$ are independent of $C$ itself. This is akin to testing for the partial correlation between valuation and price to be 0.

The results that we obtain from the data are as follows:

itemize$\Omega_3$ provides a reasonably good fit, yielding the value of the multiple $R^2$ to be $0.75$ and $0.64$ respectively for the two equations. The regression diagnostics checks were satisfied. • The correlation between the residuals of the two regression equations of $\Omega_3$ came out to be a whopping 0.859; which cannot be attributed to sheer chance with such a large sample size (a sample size of about $19,000$).

This has widespread implications. First of all, this substantiates that valuation of the auctioneers is not as harmless as providing an educated guess for the price, but rather have a causal impact on the final price. This occurs, as the base price is visible to all potential buyers. Thus, to automate the process, whatever procedure we propose cannot be held against the standard with how the valuation predicts price, as the predictability is nested as a causal impact. Hence, to automate the process, a possible change in the auction mechanism is required, and the premises of checking the success of any such alternate mechanism with this data (where valuation has had a causal impact) would be flawed.

Causal Analysis with Three-Stage Latent Hierarchical Causal Model

The linear analysis provides substantial evidence in favor of the price- valuation linkage. But this may be confounded with many other interactions present in the model. Recall that we have $6$ clusters for types of tea grades and $7$ clusters for source or tea gardens from where the packet came from. Together, considering their interactions, we have a $6\times 7 = 42$ element matrix. These cells are the hidden states in the model and can be interpreted as the proxies for demands for different types of tea in the market. These states, along with several other predictors will generate a common value about a single grade, garden, and week combination, which shall be the true value of the tea packet to both the auctioneers and the buyers. Finally, these prices and valuations will be characterized based on the common value and shall further be influenced by the particular volume of the lot, which would yield the final observations. The model that describes such a situation most closely is a variant of the Linear dynamical system model, popularly associated with Kalman filter, as described in kalman1960new and kalman1961new.

The mathematical framework is as follows;

equation[equation omitted — 125 chars of source]

where $Z_t$ is a $42\times 1$ vector, denoting the market condition. Here, no intercept term is used, since we wish the matrix $F$ to be interpreted as a transition matrix over the market condition vector $Z_t$. One can simply identify $Z_t$ to be a non-deterministic linear dynamic system. Let, $Z_t$ be denoted symbolically as;

$$ Z_t =

bmatrix[bmatrix omitted — 90 chars of source]

$$

where $a_{ij, t}$ is the latent market condition for the demand of tea in the $i$-th grade cluster and $j$-th source cluster, at time $t$. The next level of the model is;

multline[multline omitted — 203 chars of source]

where $W_{gt}$ is a scalar, which denotes the common value of the tea lot at the combination $g$ = (Garden, Grade) at time $t$, which depends on its previous observation, the current market state $Z_t$ and some exogenous control variables $X_{gt}$. In this case, the vector $G_g$ has a special structure such that;

$$G_g Z_t = \beta_1 a_{i_0, j_0, t} + \beta_2 \sum_{i \neq i_0} a_{i, j_0, t} + \beta_3 \sum_{j \neq j_0} a_{i_0, j, t}$$

where $\beta_1, \beta_2$ and $\beta_3$ are parameters to be estimated. Here, $g$ is a grade and garden combination such that the grade belongs to the $i_0$-th cluster, and the garden belongs to the $j_0$-th cluster of the source. This special structure means that the common value for a tea packet depends on the market condition of demand for that particular type of tea, as well as the market condition of its available substitutes, which shares either the same tea grade or the same source garden, as a potential cause for substitutability.

And finally, we have the model for the observations;

multline[multline omitted — 249 chars of source]

where $y_{igt}$ is the actual bivariate observation of Price and Valuation of $i$-th repeated measure in $g$-th group combination in $t$-th time, while $u_{igt}$ is some more exogenous variables, whose influences are incorporated only in the final stage and;

$$\mathbf{1}_2 =

bmatrix[bmatrix omitted — 22 chars of source]

$$

The reason we need to deviate from the standard Simple Adaptive Control Model meyntweedie (Page 40) or the popularly known Kalman Filter model (which allows for only two indices), is that we have more than one observations which are manifestations of the same state (i.e., three indices), and hence, a direct influence of the states do not account for the variability within the observations of the same state.

In the above specification of our three-stage model, we only observe the variables $y_{igt}, u_{igt}$ and $X_{gt}$, and the latent variables $Z_t, W_{gt}$ are unobservable. Hence, we can characterize the model with the specification of all the parameters and the unobservable latent variables, namely by the list of elements $(Z_t, W_{gt}, F, Q, \Phi_0, \Phi, G_g, H, R, \Gamma, S)$. Unfortunately, the above model is not identifiable, as the new set of elements given by $(\alpha Z_t, W_{gt}, F, \alpha^2 Q, \Phi_0, \Phi, \frac{1}{\alpha} G_g, H, R, \Gamma, S)$, also result in the exact same model. The crucial reason for this unidentifiablity is that the first equation contains no observable variable. For this reason, we require to pose a constraint on the model by specifying $\Vert Z_t \Vert = 1$, i.e. the vector $Z_t$'s are normalized for any $t = 0, 1, \dots T$.

commentTo perform the estimation of the parameters in the above model, we first consider the likelihood function conditional on the observed data; \begin{align*} \mathcal{L} & = f(y \mid W, u)\times f(W \mid Z, X) \times f(Z)\\ & \propto \prod_{i, g, t} {\vert S\vert}^{-1/2} \exp\left[-\frac{1}{2} (y_{igt} - W_{gt}\mathbf{1}_2 - \Gamma u_{igt} )^{\top}S^{-1} (y_{igt} - W_{gt}\mathbf{1}_2 - \Gamma u_{igt} )\right] \times \\ & \qquad \prod_{g, t} {\vert R\vert}^{-1/2} \exp\left[-\frac{1}{2} (W_{gt} - \Phi_0 - \Phi W_{g(t-1)} - G_g Z_t - HX_{gt} )^{\top}\right. \\ & \qquad \qquad \qquad \qquad \qquad \bigg. R^{-1} (W_{gt} - \Phi_0 -\Phi W_{g(t-1)} - G_g Z_t - HX_{gt} ) \bigg] \times \\ & \qquad \prod_t {\vert Q\vert}^{-1/2} \exp\left[-\frac{1}{2} (Z_{t} - FZ_{(t-1)} )^{\top}Q^{-1} (Z_{t} - F Z_{(t-1)} ) \right] \end{align*} However, this full conditional likelihood also contains the unobservable parameters, hence, the maximum likelihood estimate of parameters of the model would also depend on those unobservables. But if those unobservables are known, then the process of estimating the parameters is simply estimation of $3$ independent regression models, which can readily be estimated through means of Ordinary Least Squares (OLS). Since, the unobservables are not known, we can obtain maximum likelihood estimates of those unobservables conditional to the value of the parameters. Therefore, to obtain the estimates of the unobserved latent variables $Z_t$ and $W_{gt}$'s, we simply have to minimize the following objective function with respect to them. \begin{multline} \mathcal{F} = \sum_{t} (Z_{t} - FZ_{(t-1)} )^{\top}Q^{-1} (Z_{t} - F Z_{(t-1)} ) + \sum_{g, t} \frac{1}{R}(W_{gt} - \Phi_0 - \Phi W_{g(t-1)} - G_g Z_t - HX_{gt} )^2\\ + \sum_{i, g, t} (y_{igt} - W_{gt}\mathbf{1}_2 - \Gamma u_igt )^{\top}S^{-1} (y_{igt} - W_{gt}\mathbf{1}_2 - \Gamma u_igt ) \end{multline} Now, we first differentiate above with respect to $Z_t$; \begin{multline} \frac{\partial{\mathcal{F}}}{\partial Z_t} = 2 Q^{-1} (Z_t - F Z_{(t-1)}) + 2 F^{\top} Q^{-1} (Z_{(t+1)} - F Z_t)\\ + \frac{2}{R} \sum_g G_g^{\top} (W_{gt} - \Phi_0 - \Phi W_{g(t-1)} - G_g Z_t - HX_{gt} ) \end{multline} and differentiating with respect to $W_{gt}$; \begin{multline} \frac{\partial{\mathcal{F}}}{\partial W_{gt}} = \frac{2}{R} (W_{gt} - \Phi_0 - \Phi W_{g(t-1)} - G_g Z_t - HX_{gt} ) \\ + \frac{2}{R} \Phi(W_{g(t+1)} - \Phi_0 - \Phi W_{gt} - G_g Z_{(t+1)} - HX_{g(t+1)} ) + \sum_i \mathbf{1}_2^{\top} S^{-1} (y_{igt} - W_{gt}\mathbf{1}_2 - \Gamma u_{igt}) \end{multline} Finally, we have the following estimating equations, which are obtained by setting the above derivatives equal to $0$: \begin{multline} Z_t = \left[Q^{-1} - F^{\top}Q^{-1}F + \frac{1}{R} \sum_{g}G_g^{\top}G_g\right]^{-1} \Bigg[ Q^{-1}FZ_{(t-1)} - F^{\top}Q^{-1}Z_{(t+1)}\\ \left. - \frac{1}{R} \sum_{g}G_g^{\top} (W_{gt} - \Phi_0 - \Phi W_{g(t-1)} - HX_{gt}) \right] \end{multline} and \begin{multline} W_{gt} = \dfrac{1}{\left(\dfrac{1 + \Phi^2}{R} - r_{gt}\mathbf{1}_2^{\top} S^{-1}\mathbf{1}_2 \right)} \times \left[ \dfrac{\Phi_0 (1+\Phi)}{R} + \dfrac{\Phi (W_{g(t-1)} - W_{g(t+1)})}{R} + \dfrac{G_g (Z_t + \Phi Z_{(t+1)})}{R} \right.\\ \left. + \dfrac{H (X_{gt} + \Phi X_{g(t+1)}) }{R} - \mathbf{1}_2^{\top} S^{-1}\left(y_{igt} - \Gamma u_{igt}\right) \right] \end{multline} Since, equation (ref) contains $W_{gt}$ in its right hand side, and equation (ref) contains $Z_t$ in its right hand side, we should try to solve these simultaneous linear equations for both $Z_t$ and $W_{gt}$ together. Putting, the value of $W_{gt}$ as obtained from equation (ref) into equation (ref), we obtain; \begin{multline*} \left[Q^{-1} - F^{\top}Q^{-1}F + \frac{1}{R} \sum_{g}G_g^{\top}G_g\right] Z_t\\ = \Bigg[ Q^{-1}FZ_{(t-1)} - F^{\top}Q^{-1}Z_{(t+1)} - \frac{1}{R} \sum_{g}G_g^{\top} \Bigg\{ \left(\dfrac{1 + \Phi^2}{R} - \dfrac{r_{gt}}{S}\right)^{-1} \times \Bigg( \dfrac{\Phi_0 (1+\Phi)}{R} + \dfrac{G_g (Z_t + \Phi Z_{(t+1)})}{R} \\ + \dfrac{\Phi (W_{g(t-1)}- W_{g(t+1)})}{R} + \dfrac{H (X_{gt} + \Phi X_{g(t+1)}) }{R} - \dfrac{y_{igt} - \Gamma u_{igt}}{S} \Bigg) - \Phi_0 - \Phi W_{g(t-1)} - HX_{gt}\Bigg\} \Bigg] \end{multline*} Now, this equation can be simplified to obtain the updated value of $Z_t$ given by; \begin{align} \begin{split} Z_t & = \left[Q^{-1} - F^{\top}Q^{-1}F + \frac{1}{R} \sum_{g}G_g^{\top}G_g + \frac{1}{R^2} \sum_{g}G_g^{\top} \left(\dfrac{1 + \Phi^2}{R} - r_{gt} \mathbf{1}_2^{\top} S^{-1} \mathbf{1}_2 \right)^{-1} G_g \right]^{-1}\\ & \qquad \Bigg[ Q^{-1}FZ_{(t-1)} - F^{\top}Q^{-1}Z_{(t+1)} - \frac{1}{R} \sum_{g}G_g^{\top} \Bigg\{ \left(\dfrac{1 + \Phi^2}{R} - r_{gt} \mathbf{1}_2^{\top} S^{-1} \mathbf{1}_2 \right)^{-1} \times\\ & \qquad \Bigg( \dfrac{\Phi_0 (1+\Phi)}{R} + \dfrac{G_g (\Phi Z_{(t+1)})}{R} + \dfrac{\Phi (W_{g(t-1)}- W_{g(t+1)})}{R} + \dfrac{H (X_{gt} + \Phi X_{g(t+1)}) }{R}\\ & \qquad - \mathbf{1}_2^{\top} S^{-1} \left(y_{igt} - \Gamma u_{igt} \right) \Bigg) - \Phi_0 - \Phi W_{g(t-1)} - HX_{gt}\Bigg\} \Bigg] \end{split} \end{align} Therefore, using equation (ref) one firstly computes the updated value of $Z_t$, then due to the constraint $\Vert Z_t \Vert = 1$, we normalize it and then equation (ref) can be used to compute the updated value of $W_{gt}$ using the updated value of $Z_t$ in its right hand side. The complete estimation algorithm in given in Algorithm (ref). \begin{algorithm}[ht] \SetAlgoLined \SetKwInOut{Input}{Input} \SetKwInOut{Output}{Output} \Input{$X_{gt}, u_{igt}, y_{igt}$} \Output{The estimates $F, Q, \Phi_0, \Phi, \beta_1, \beta_2, \beta_3, H, R, \Gamma, S$} Set some values of $Z_t$ and $W_{gt}$ for each $g = 1, 2, \dots N_t; t = 1, 2, \dots T$\; \Repeat{convergence}{ $F \leftarrow \left[\sum_{t} Z_{t} Z_{(t-1)}^{\top} \right] \left[\sum_{t} Z_{(t-1)} Z_{(t-1)}^{\top} \right]^{-1}$\; $Q \leftarrow \dfrac{1}{T-1} (Z_t - FZ_{(t-1)})(Z_t - FZ_{(t-1)})^{\top}$\; Estimate $\Phi_0, \Phi, \beta_1, \beta_2, \beta_3, H$ from equation (ref)\; $R \leftarrow \dfrac{1}{\sum_t N_t - p} (W_{gt} - \Phi_0 - \Phi W_{g(t-1)} - G_g Z_t - HX_{gt})(W_{gt} - \Phi_0 - \Phi W_{g(t-1)} - G_g Z_t - HX_{gt})^{\top} $, where $p$ is the total number of parameters in equation (ref)\; $\Gamma \leftarrow \left[\sum_{t} (y_{igt} - W_{gt}\mathbf{1}_2 ) u_{igt}^{\top} \right] \left[\sum_{t} u_{igt} u_{igt}^{\top} \right]^{-1}$\; $S \leftarrow \dfrac{1}{\sum_{g,t} r_{gt}} (y_{igt} - W_{gt}\mathbf{1}_2 - \Gamma u_{igt})(y_{igt} - W_{gt}\mathbf{1}_2 - \Gamma u_{igt})^{\top}$\; \For{$t = 1, 2, \dots T$}{ Use equation (ref) to update the value of $Z_t$\; $Z_t \leftarrow Z_t / \Vert Z_t \Vert$\; \For{$g = 1, 2, \dots N_t$}{ Use equation (ref) to update the value of $W_{gt}$\; } } } \caption{Algorithm for estimation of Three Stage Latent Hierarchical (TSLH) Model} \end{algorithm}

Figure (ref) shows the Causal DAG diagram for the above three-stage latent hierarchical (TSLH) model, for a fixed time point $t$, the defining SCM for this are given by equations (ref), (ref) and (ref). Note that, there are $4$ exogenous variables at the second stage (as noted from the four parameters $H_1, H_2, H_3$ and $H_4$) and only one exogenous variable in the third stage. We use Gibb's sampler to obtain the estimates, as discussed in the following subsection. The description of these parameters, along with the estimated value from the dataset is given in table (ref).

figure[figure omitted — 1,889 chars of source]

Gibbs Sampling Conditional Derivations

We begin by writing the likelihood for the Three Stage Latent Hierarchical (TSLH) model, upto a proportionality constant.

align*[align* omitted — 639 chars of source]

Before obtaining the individual conditional distributions, the following observation will come in handy:

Note that, when $Y\sim \mathcal{N}_k(\mu,\Sigma)$ as in regression model, then $$f(\mathbf{y})\propto \exp\left(-\dfrac 12Y^\top \Sigma^{-1} Y + Y^\top \Sigma^{-1}\mu \right)$$ Therefore, we can simply identify the normal distribution based on these coefficients $\Sigma^{-1}$ and $\Sigma^{-1}\mu$, which is a reparametrization of the parameters of normal distribution. This reparametrization shall be helpful in identifying the conditional distributions in the subsequent calculations.

We first obtain the conditional distributions for the latent variables, $Z_t$ and $W_{gt}$ respectively.

multline*[multline* omitted — 378 chars of source]

which is a normal distribution with the above reparametrization where, $\Sigma^{-1}$ is the coefficient of quadratic term, and $\Sigma^{-1}\mu$ is the coefficient of the linear term.

Now, with $W_{gt}$, we have;

multline*[multline* omitted — 490 chars of source]

This leads to another normal distribution. Next, for the parameters,

$$Q \mid \text{rest} \propto \vert Q\vert^{-T/2} \exp\left[ -\dfrac{1}{2} \text{tr}\left( \left(\sum_t \xi_t \xi_t^{\top} \right) Q^{-1} \right) \right]$$

which leads to $\mathcal{W}^{-1}\left[ \sum_t \xi_t \xi_t^{\top}; T-43 \right]$ distribution, where $\mathcal{W}^{-1}$ is used to denoted Inverse Wishart distribution.

Similarly,

$$R \mid \text{rest} \propto R^{-\sum_{t}N_t /2} \exp\left[ -\dfrac{1}{2} \left(\sum_{g,t} e_{gt} e_{gt}^{\top} \right) R^{-1} \right]$$

which is same as the density function of Inverse Gamma distribution with shape $\alpha = \sum_t \dfrac{N_t}{2} - 1$, scale parameter $\beta = \dfrac{1}{2} \left(\sum_{g,t} e_{gt} e_{gt}^{\top} \right)$, upto a proportionality constant.

And finally,

$$S \mid \text{rest} \propto \vert S\vert^{-N/2} \exp\left[ -\dfrac{1}{2} \text{tr}\left( \left(\sum_{i, g, t} \epsilon_{igt} \epsilon_{igt}^{\top} \right) S^{-1} \right) \right]$$

which again leads to $\mathcal{W}^{-1}\left[ \sum_{i, g, t} \epsilon_{igt} \epsilon_{igt}^{\top}; N-3 \right]$ distribution.

Continuing,

align*[align* omitted — 473 chars of source]

therefore,

$$F \mid \text{rest} \sim \mathcal{MN}_{42\times 42}\left( \left[ \sum_t Z_{t-1}Z_{t-1}^{\top}\right]^{-1} \left[ \sum_t Z_t Z_{t-1}^{\top}\right], I, \left[ \sum_t Z_{t-1}Z_{t-1}^{\top}\right]^{-1} Q \right)$$

where $\mathcal{MN}$ stands for the matrix normal distribution. In other words, since the covariance matrix between the rows of the $F$ is $I$, the identity matrix, hence we can generate the rows of $F$ independently from multivariate normal distributions with mean vectors same as the rows of the mean matrix, and the same covariance matrix $\left[ \sum_t Z_{t-1}Z_{t-1}^{\top}\right]^{-1} Q$.

On a similar note,

align*[align* omitted — 623 chars of source]

Therefore,

$$\Gamma \mid \text{rest} \sim \mathcal{MN}_{2\times 2}\left( \left[\sum_{i, g, t} u_{igt} u_{igt}^{\top}\right]^{-1} \left[\sum_{i, g, t} \left( y_{igt} - W_{gt}\mathbf{1}_2 \right)u_{igt}^{\top}\right], I, \left[\sum_{i, g, t} u_{igt} u_{igt}^{\top}\right]^{-1}S \right)$$

To get conditional distribution of parameters corresponding to stage 2 of the model, we assume that, $G_g Z_t = \beta_1 \tilde{Z}_{1t} + \beta_2 \tilde{Z}_{2t} + \beta_3 \tilde{Z}_{3t}$, with $\tilde{Z}$ being the proper linear combination of latent state $Z_t$ that affects the common value $W_{gt}$. Let us also denote the vector of parameters,

$$\theta =

bmatrix[bmatrix omitted — 66 chars of source]

$$

Based on this, we have;

$$\theta \mid \text{rest} \sim \mathcal{MVN}\left( A^{-1}b; A^{-1}R \right)\text{ where }b = \left( \sum_{g, t} W_{gt}, \sum_{g, t} W_{gt} W_{g(t-1)}, \sum_{g, t} W_{gt} \tilde{Z}_{1t}, \dots \right)^\top $$

and

$$A =

bmatrix[bmatrix omitted — 453 chars of source]

$$

where $\mathcal{MVN}$ denotes the multivariate normal distribution.

Estimated values and implications

The performance of the estimated model has been shown in Figure (ref) and in Figure (ref). Some of the residual diagnostics are shown in Figure (ref) and Figure (ref), from which it is obvious that the residuals in logarithm scale follow an approximate normal distribution, other than some outlying values in the tail, as well as the residuals in the original scale of price, shows a histogram of leptokurtic distribution which closely resembles a lognormal one. The estimated values of the parameters of TSLH model are shown in the table (ref).

longtable[longtable omitted — 4,357 chars of source]
figure[figure omitted — 587 chars of source]
figure[figure omitted — 585 chars of source]
figure[figure omitted — 570 chars of source]
figure[figure omitted — 549 chars of source]

To emphasize the goodness of fit for the Three Stage Latent Hierarchical (TSLH) model we obtain the following:

enumerate• The residuals in valuation predicting component of the TSLH model makes approximately $38$ rupees of error with $80\%$ confidence and about $59$ rupees of error with $95\%$ of confidence. • The residuals in price predicting component of the TSLH model makes approximately $22$ rupees of error with $80\%$ confidence and about $29$ rupees of error with $95\%$ of confidence.

The results we obtain have the following notable implications:

enumerate• Since $\beta_1$ is positive, it appears that the more is the market demand for a particular type of tea, the more is its value to the buyers. However, its small value (inclusion of 0 in confidence interval) suggests that this dependence with the underlying market condition is nonetheless small. • Since $\beta_2$ and $\beta_3$ are estimated to be negative, it appears that the substitution effect is present in this scenario. The more is the market demand for the substitutes of a tea type, the less is its utility to the buyers, who are going to ultimately sell it to the consumers. • Negative values of $H_1, H_2$, and $H_3$ show that, if one type of tea has already dispatched in the market by a large volume in the recent past, then their utility to the buyers of the auction is going to be less. • An interesting thing to note is the opposite behavior of $\Gamma_p$ and $\Gamma_v$ when we allow an unrestricted $\Gamma_{pv}$ in the model. Generally, we expect that larger volumes of tea packets to be sold at smaller per unit price, which is in accordance with a negative $\Gamma_p$ value. However, if the valuation is indeed nothing but an educated guess for the price, then the behavior for $\Gamma_p$ and $\Gamma_v$ should be identical. The fact that $\Gamma_v$ is positive suggests that the auctioneer tries to balance for the act of the decrease in per-unit price with larger packet volumes, by setting its valuation at a price larger than the market would expect it to be. • A large positive value of $\Gamma_{pv}$ suggests a very strong positive dependence on the prices of tea packets on the auctioneers' valuation of the tea packets. However, this dependence contains both the direct dependence of price on the valuation and an indirect dependence of price through the mediating effect of the spurious latent variable $W_{gt}$.

Test of causality

Now, to answer the question about whether there is a significant direct effect from valuation to the price, (i.e. whether the dotted arrow in Figure (ref) exists or not) a very general approach is the conditional independence test, as discussed in pearl2016causal and AnIntroductiontoCausalInference. As shown in Figure (ref), $W_{gt}$ and $u_{igt}$ creates a fork with $\log(\text{Valuation}_{igt})$ and $\log(\text{Price}_{igt})$ nodes in the DAG. However we shall require the following theorem:

theoremThe estimated $\Gamma_{pv}$ based on equation (ref) where $W_{gt}$ is substituted by the samples obtained from the Gibbs sampler, converges in distribution to the $\Gamma_{pv}$ estimated using the true value of the latent variables (although unknown).
proofWe shall use $f(\cdot)$ and $F(\cdot)$ as a generic symbol for density and distribution function of a probability distribution, respectively. Let, $W_{gt, k}^{(r)}$ be the samples taken from the simulated $k$-th Markov chain in Gibbs sampling, after $r$ iterations (i.e. after $r$ transition steps of the Markov chain), for the latent variable $W_{gt}$. It is well known that under certain regularity conditions such as positivity and connected of the conditional densities the simulated Markov chain will converge to a stationary distribution which is the desired posterior distribution, as shown in robert2013monte. Therefore, as $r\rightarrow \infty$, we should have $W_{gt, k}^{(r)} \xrightarrow{\mathcal{L}} W_{gt, k}$, where $W_{gt, k}$ is a random variable following the posterior distribution with density $f(W_{gt} \mid \text{data})$. Here, the notation $\xrightarrow{\mathcal{L}}$ is used to denote convergence in distribution for random variables. Hence, the joint distribution of the observed data together with the posterior samples from Gibbs sampling, converges to the joint distribution of the observed data together with unobserved latent variables. \begin{align*} & F\left(y_{igt}, u_{igt}, W_{gt, k}^{(r)}; \forall i, g, t\right) \\ = \quad & F\left(y_{igt}, u_{igt}; \forall i, g, t\right) F\left(W_{gt}^{(r)} \mid y_{igt}, u_{igt} ; \forall i, g, t\right)\\ \xrightarrow{\mathcal{L}} \quad & F\left(y_{igt}, u_{igt}; \forall i, g, t\right) F\left(W_{gt} \mid y_{igt}, u_{igt} ; \forall i, g, t\right)\\ = \quad & F\left(y_{igt}, u_{igt}, W_{gt, k}; \forall i, g, t\right) \end{align*} where we use the fact that $W_{gt}^{(r)} \xrightarrow{\mathcal{L}} W_{gt}$ means the distribution function converges pointwise. Now the proof follows once we note that, $\Gamma_{pv}$, as a regression estimate, is a continuous function of $y_{igt}, u_{igt}$ and $W_{gt}$. Therefore, by continuous mapping theorem, \begin{equation} \left[\sum_{i, g, t} u_{igt} u_{igt}^{\top}\right]^{-1} \left[\sum_{i, g, t} \left( y_{igt} - W_{gt, k}^{(r)}\mathbf{1}_2 \right)u_{igt}^{\top}\right] \xrightarrow{\mathcal{L}} \left[\sum_{i, g, t} u_{igt} u_{igt}^{\top}\right]^{-1} \left[\sum_{i, g, t} \left( y_{igt} - W_{gt}\mathbf{1}_2 \right)u_{igt}^{\top}\right] \end{equation} Let, $\hat{\Gamma}$ be the regression estimate if the true value of $W_{gt}$ were known, while $\tilde{\Gamma}_{k}^{(r)}$ be the regression estimate where $W_{gt}$ is substituted by the posterior sample $W_{gt, k}^{(r)}$. Then, equation (ref) simply means $\tilde{\Gamma}_{k}^{(r)} \xrightarrow{\mathcal{L}} \hat{\Gamma}$, and hence particularly one entry of the matrix $\tilde{\Gamma}_{pv, k}^{(r)} \xrightarrow{\mathcal{L}} \hat{\Gamma}_{pv}$. Thus, the result follows.

There are particularly two remarks to be made relating to the above theorem.

enumerate$\tilde{\Gamma}_{pv, k}^{(r)}$ is not the posterior samples obtained from the Gibbs sampling, but rather a regression estimate based on the posterior sample $W_{gt, k}^{(r)}$. • Using this theorem, we know that under the null hypothesis that valuation has no direct causal effect on price in the TSLH model, a test for $H_0: \Gamma_{pv} = 0$ can be conducted using the regression estimate based on the proxies of the latent variables obtained from Gibb's sampler values. In other words, for each Markov chain, we can obtain samples $\tilde{\Gamma}_{pv, k}^{(r)}$, and as $r \rightarrow \infty$, these samples behave like independent and identically distributed samples from the distribution of $\hat{\Gamma}_{pv}$. Therefore, we can simply use the Hybrid confidence sets obtained from the Gibbs sampling to test our hypothesis.

To test the causality from $\log(\text{Valuation}_{igt})$ to $\log(\text{Price}_{igt})$, we shall require a conditional independence test between these two variables conditioned on the value of $W_{gt}$ and $u_{igt}$. In the given model, such independence would hold if and only if $\Gamma_{pv} = 0$. Therefore, in view of the above theorem, using the posterior samples $W_{gt, k}^{(r)}$, we obtain the estimate of $\hat{\Gamma}$ as $1.006317$, and the $95\%$ confidence set turns out to be $(0.981592, 1.020155)$, which does not contain $0$, thereby, showing sufficient evidence against the null hypothesis of conditional independence.

On the other hand, based on the fitted model, let us consider the residuals from Valuation predicting component, and the residuals from Price predicting component (leaving Valuation as an explanatory variable), and denote their product moment correlation as $r_{pv}$. In other words,

$$r_{pv} = \text{cor}\left( \log(\text{Valuation}_{igt}) - W_{gt} - \Gamma_v \log(\text{Volume}_{igt}), \log(\text{Price}_{igt}) - W_{gt} - \Gamma_p \log(\text{Volume}_{igt}) \right)$$

Tracking this correlation where the latent variables and parameters are substituted by posterior samples obtained from the Gibbs sampler, we obtain the posterior mean of $r_{pv}$ as $0.7819911$. In contrast to that, the Pearson's correlation coefficient between the logarithms of those variables valuation and price is $0.9573915$, and the Pearson's correlation coefficient between those variables valuation and price without any transformation is $0.9609977$. Note that, the correlation between the residuals obtained from the causal linear model described before was $0.859$. Therefore, we find that, the linear model was enough to structurally model some of the dependence, while the TSLH model was further able to reduce the correlation by explaining temporal dependence structure within the data. However, there was still unexplained correlation, which was simply a manifestation of the causal relationship between auctioneers' valuation and the ultimate selling price at the auction.

Remarks on Automation of the Process and Conclusion

The preceding sections show us the utmost significance of the manual valuation of the tea packets that come in, in predicting the final price level. Thus the hope of automating the entire procedure seems unrealistic.

However, since this valuation is based on the inherent characteristics of the tea dust packets, and the volume of packets that arrive, we strongly believe that the practice of using valuations to set base prices can be done away with. George and Hui, in their paper optimal, provide an ingenious way to estimate demand in the auction market, under the Independent Private Value Model second-price auctions. A generalization of this method, in this regard, to the Common Value (CV) auction_book case, where the optimal symmetric bidding strategies are not the bidder signals themselves, but a monotonic function of their signals, may be helpful in our case. Then, knowing the distribution of the bidder values, an optimal reserve price may be set to maximize the ex-ante expected revenue, which is a function of this distribution. Levin and Smith, in their paper disproof, have shown that under the non-IPV case, the optimal reserve price for the seller converges to her true value - here it's her manufacturing costs. Hence if the pool of bidders grow, then it would be safe for the seller to set the reserve price at her manufacturing costs. Collusion among the bidders is often a very practical problem to ponder about, and most methodologies fail under scenarios not robust to such behavior. For example, in second-price auctions, one source of asymmetric equilibrium, (when the distribution has a support $[0,\omega]$) is for one bidder to bid $\omega$ and the others to bid $0$(or the minimum possible price). This is often realized in real-life scenarios, e.g. spectrum auctions. This is a possible solution here, and given that several auctions occur regularly in this market, bidders can sequentially alternate the role of the highest bidder, and thus can all be better off, at the cost of the seller. Our pricing model provides a way to detect such behavior on the part of the bidders. Since our pricing model with valuations provides an excellent fit to the true prices, this can be used to detect collusion. As in the case of collusion, the final price of the transaction would be low compared to the expected transaction price, a large deviation from the predictions would indicate the presence of such collusion. Our prices, with the estimated parameters, approximately follow a normal distribution, thus a low p-value from this distribution could be used as an indication of collusion of bidders. Further research on the aforesaid aspects could bring exciting breakthroughs in the path of automation.