1. Introduction
India is primarily an agriculturebased country and its economy largely depends upon agriculture. According to World Bank, approximately 60 percent of India’s land area is used for agricultural purpose making India the second largest in terms of agricultural land availability. Major chunk of this agricultural land is rainfed (60% of total agricultural land) and India ranks first among the rainfed agricultural countries. Due to unpredictable monsoons there is constant threat to sustainable agricultural production. Presently, contribution of agriculture is about one third of the national GDP and provides employment to over 80% of Indian population in agriculture and allied activities. Consequently, sustainable development of India largely depends upon the development of agriculture. The agricultural production information is very important for planning and allocation of resources to different sectors of agriculture. This is more important for Karnataka, the second largest droughtprone state in India after Rajasthan. The state of Karnataka largely depends on weather conditions for effective farming with more than 75 percent of its arable land in the rainfed (farming by rain water only) regions. Thus, minimizing the impact of natural disasterrelated crop losses, particularly from drought, is therefore a major public policy objective for Karnataka government.
The Karnataka Agriculture Price Commission (KAPC) was set up in 2013 to address several agrirelated issues. It works to ensure a) maximum share of consumer price, b) achieve sustainable development in the field of agriculture in the State, c) provide remunerative price for farm produce and d) provide a suitable marketing system^{1}^{1}1http://kapricom.org/downloads/about/termsOfRefen.pdf.
Red gram is commonly known as Tur or Arhar (Pigeon pea) in India and is the second most important crop in the pulses category in the country after gram (chana). The ability of red gram to produce high economic yields under soil moisture deficit makes it an important crop in rainfed and dryland agriculture. India contributes for nearly 65% of world’s total red gram production. However, it is gaining importance in African countries due to its adaptability to limited moisture conditions and also favorable nitrogen fixation properties exhibited by the crop.
The government of India publishes minimum support prices (MSP) for important agriculture commodities, and Tur is an important crop in that respect. For the last few years Tur has seen huge market disruptions owing to poor prediction of arrivals as well as corresponding MSP. There has been sharp variations in prices on a year to year basis and it has led to lot of focus in the prediction of the arrivals to market, both from the government side and also from the agroresearch community perspective. The government of Karnataka along with KAPC has studied the commodity arrival patterns and the subsequent problem of price prediction problem for several years. Arrival prediction helps in identifying the predicted yield for a particular commodity, which further helps in planning the price related interventions. KAPC is collaborating with Microsoft to understand and predict these patterns. The model and the subsequent insights described in this work were presented to the state government at the end of last year.^{2}^{2}2https://www.indiatoday.in/ptifeed/story/karnagovtsignsmouwithmicrosoftfortechnologyoriented107403920171027
https://www.techgig.com/technews/microsoftaitohelpkarnatakafarmersgethighercropyields137863
http://www.thehindubusinessline.com/economy/farmerslooktoharvestthefruitsofai/article9928335.ece
http://indianexpress.com/article/india/karnatakagovtinksmouwithmicrosofttouseartificialintelligencefordigitalagriculture4909470/
In this work we describe a market wise commodity arrivals prediction (for total arrivals in a month to a market) over several primary markets and also at the entire state level. We present results against standard methods and show that the remote sensing based method performs better than just market based timeseries prediction models or other popular machine learning (ML) based methods.
1.1. Related Work
Economists have historically dealt with small sets of data, but this constraint is rapidly changing in the current day and age as new and more detailed data becomes easily available. With the advent of remote sensing, the amount of data that economists can look at went to several gigabytes and ranging over millions of observations. One can now analyze large chunks of data with very little biases / corrections / errors introduced by human manipulations. According to (Donaldson2016, )
, the main advantages of remote sensing data to economists can be classified into different categories: 1) access to information which is difficult to obtain by other means; 2) unusually high spatial resolution.
The first advantage arises from the fact that remote sensing satellites can collect panel data at low marginal cost, repeatedly over large time epochs, and usually at very large scale. This data can then be treated as a proxy for a wide range of difficult to measure characteristics. Many topics of potential interest have already been measured remotely and used in other fields such as urban development, building types, road characteristics, pollution indices, beach quality and many more. The ground data needed to study these fields in depth would be prohibitively expensive, and difficult to measure with desirable accuracy. Additionally, there are scenarios when the official government counterparts of some remotely sensed statistics (such as pollution or forestry) may be subject to human manipulation.
The second advantage of remote sensing data sources is that they are typically available at a substantially higher degree of spatial resolution than are traditional data. Much of the publicly available satellite imagery used by economists provides readings for each of the hundreds of billions of 30m30m grid cells of land surface on Earth. The freely available MODIS data (MODIS, ; MODIS2, ) gives resolutions as 250m, 500m and 1000m. These resolutions are mostly sufficient for econometric studies. Even higher resolution data can be obtained for very targeted applications / studies.
Given the rate at which algorithms for crop classification and yield measurement have improved in recent years, future applications of satellite data are likely to be particularly rich in the agricultural arena. Different bands, and combinations thereof, have various useful properties. Plants reflect different sets of frequencies, at different rates, for different stages of their life cycle. For this reason, functions of reflectance in specific portions of the visible and infrared spectra can provide information about vegetation growth. This insight is used to produce the commonly used Normalized Difference Vegetation Index (NDVI) (Diorio1989, ). Infrared data can also be used to measure temperature and, indirectly, precipitation when applied to clouds (for example, (Novella2013, )). Different techniques have been proposed to merge remote sensing data, more specifically NDVI with yield and crop parameter prediction. Time series based methods for yield prediction were proposed by Rembold et al. (Felix2013, ). Supervised classification techniques were proposed by Doraiswamy et al. (Doraiswamy04, )
. Linear regression based crop coefficient prediction models were proposed by Kamble et al.
(Kamble2013, ). In a recent paper Aviv and LundsgaardNielsen (Aviv2017, ) propose a method to predict soy yield by coupling remote sensing data with soil measurement data for 58 locations.You et al. (You2017, ) propose a new method where they go even beyond NDVI and work directly with the band images obtained from satellites to estimate cop yields. Their method is shown to work good for countries like USA due to ready availability of large amounts of clean data. We concentrate solely on data from the state of Karnataka (KA), India, which is available from 2013 onwards on a monthly basis. The amount of data is not enough to train the data intensive neural method described in (You2017, ). Yield estimation is a difficult problem which needs clean data both for modeling as well as validation purposes. You et al. (You2017, ) use USDA published data as the ground truth for validating their models. Similar quality data is hard to obtain for developing countries. Moreover, the farm sizes as well as cropping patterns are very different in the developing countries as compared to countries like USA.
In this paper, our contribution is proposed framework to work with limted timeframe data at large scale, a situation often faced in delevoping countries. To our knowledge no work has been done in the mentioned context, especially for any SouthAsian countries to predict market arrival and associated insights.
2. Data
2.1. Normalized Difference Vegetation Index (NDVI)
NDVI is calculated from the differential rates of visible and nearinfrared light reflected and absorbed by vegetation. Healthy vegetation absorbs most of the visible light that hits it, and reflects a large portion of the nearinfrared light. Unhealthy or sparse vegetation reflects more visible light and less nearinfrared light. Nearly all satellite based vegetation indexes employ this difference to quantify the density of plant growth on the Earth. One of the most common indexes, popular with the research community, is the Normalized Difference Vegetation Index (NDVI). NDVI is computed as the ratio between the difference and the summation of the nearinfrared (NIR) and the visible radiation (VIS) as shown in Eq. 1.
(1) 
where is the reflection for Near Infra Red bands and is the reflection for Visible spectrum. Calculations of NDVI for a given pixel (location on Earth), always results in a number that ranges from minus one (1) to plus one (+1). However, no green vegetation gives a value close to zero. A zero means no vegetation and close to +1 (0.8  0.9) indicates the highest possible density of green vegetation.
In our work, MODIS data was used for the compuatation of NDVI. MODIS uses NIR and red wavelength bands for the computation of NDVI. The spatial resolution of the data is . NDVI data is extremely noisy due to the inherent difficulty in capturing the data. We smooth the NDVI data by averaging over spatial neighborhoods as shown below.
(2) 
where is the raw NDVI value obtained from satellite images at location and time , is the neighborhood of , is the smoothed transformation output, and is the transformation function. In the current work we define as the identity function. Putting nonlinear combinations for the raw NDVI values is left as a future investigation.
We perform a sampling of the locations to further reduce redundancies in the NDVI data. We divide the geographical region of interest, state of Karnataka, into uniform blocks, and select the centroid of the blocks as the representative locations. These representative locations are selected such that they coincide with the major markets of Tur in Karnataka. Once the locations are identified we perform the transformation and smoothing of the centroid as mentioned in Eq. 2.
Commodity  Types  Commodity  Types 

Cereals  15  Drugs & Narcotics  3 
Dry Fruits  2  Fibre Crops  3 
Forest Products  8  Fruits  20 
Live Stock/Poultry  11  Oil Seeds  18 
Other  12  Pulses  21 
Spices  10  Vegetables  47 
2.2. Market Data
Commodity market data for the state of Karnataka is updated daily at the state market portal (KrishiMVahini, ). The state government publishes daily arrivals and prices data for the list of commodities as shown in Table. 1. The data available from Jan’2013 onward is reliable and relatively noise free. For each of its 82 markets, the government of Karnataka publishes the commodity variety, grade, arrivals, minimum, maximum and modal prices for each day. Out of the 170 commodities tracked by the state government, the central government of India provides minimum support price (MSP) for only 14 crops with a total of 17 varieties (MSP, ). Tur is one of the crops for which the MSP is provided and Karnataka is the second largest producer of Tur in India accounting for around 17% of total production, which is almost 11% of the Tur production across the globe. The arrivals distribution over time for some of the Tur producing regions of Karnataka are shown in Fig. 2.
3. Model
Due to historically heavy production of Tur in the northern parts of the state, this region is called the TUR bowl of Karnataka. The primary thesis of this work is the assumption that the NDVI values computed for these regions, aggregated over a few years would be a good indicator to predict the arrival of the commodity to the local markets. The change of NDVI over a calendar year for the major crop area of Kalaburgi in Karnataka is shown in Fig. 3. The green dots show higher NDVI owing to increased greenery. Yellow dots show intermediate vegetation growth over the same area and red dots show low growth at different times of the year. This pattern is observed for the Kharif season (June / July rainy season sowing). We have collected NDVI data for more than 6 Million locations distributed all over Karnataka. This data is smoothed according to Eq. 2, where the transformation function is the identity function and the free parameter is the neighborhood size . This reduces the number of NDVI data points for the entire state to Million locations. We assume that locations situated nearer to the markets provide more useful information than those situated further away. Hence, a set of locations filtered based on their proximity from the markets are selected for further processing.
Let us assume, a market at location and consider all locations such that the distance . Our method collects the NDVI values for all the locations which are close to the market (), and uses them to learn a model to predict the arrival of Tur to the market . The number of distinct locations is of the order of
K. In data cleaning exercise, we used linear inetrpolation to fill few missing observations in market data and removed high variance locatons(more than 75percentile) among the available set of locations. It was empirically observed that these high variance locations adversely effect future prediction accuracy. After data cleaning, we applied systematic sampling such that for a given market
K locations are available.Traditionally, researchers have been using principle component to reduce dimension or ridge regression to avoid noninvertible matrix, which arises in high dimension data or where predictors are likely to have high correlation among each other. In current setup, we have extremely large dimensions as a set of distinct locations, and nearby locations are likely to have highly correlated NDVI values. We evaluated the performance of these two methods and shown comparative results at the end of the section. We also considered estimating seasonality component but later discarded that thought based on arguments provided by agriculture scientists. The main concern was the fact that crop maturity depends heavily on water availability (soil moisture) and appropriate weather conditions, which does not strictly follow seasonal pattern, especially in rainfed regions where significant variations was observed in weather pattern on available last three/four years data. Limited historical data (of a few years) along with disrupted weather pattern would adversely effect any estimate of seasonal components and would only add noise in the predictions.
Based on the assumptions mentioned till this point, we explain our proposed regression framework, where first important variables are selected using penalized regression and then principle component regression is applied on selected variables to predict response variables.
Let us assume
is the vector of NDVIs for all locations
falling within a distance of market at time . The arrival of Tur for the same market at time is then predicted as an elasticnet (Zou05, ; Friedman09regularizationpaths, ) based regression outcome. We collect data for multiple time points into one single optimization and then minimize the following functional:where the input matrix is the matrix obtained by appending the NDVI data for all the locations near a market for several time instances , and the additional column is all ones to account for the bias factor . The number of rows denotes the time window in the past which is used for the regression. The second term in Eq. 3 is the penalty which penalizes nonsmooth combination vector thereby bringing correlated factors together (isl01, ). The third term with the penalty enforces sparsity in the combination vector, such that only important locations survive in the final outcome prediction. The two constraints together form the elasticnet penalty. The Jacobian for Eq. 3 with respect to the parameters can now be written as
(4) 
where gets the element wise sign of the vector argument. The update for the elements of can now be written as shown in (Friedman07pathwisecoordinate, ) as
(5) 
where
(6)  
(7) 
is the fitted value excluding the contribution from , and hence is the partial residual for fitting . is commonly called the softthresholding operator. The bias term can be obtained by solving the closed form equation
(8) 
In our formulation more weight is given to penalty , to encourage keeping as many important factors as possible. High values of
tend to remove too many locations and consequentially the derived features from these sparse locations become unstable. The instability is probably due to fact that when
is large, very few locations are picked to represent a region and therefore derived feature would be based on small sample and hence unstable.The output of the elasticnet penalty is several orders lesser in dimension than the initial size of the neighborhood . But the number of dimensions to work with is still upto a 1000 or so. We further reduce the dimension by first generating new factors as the principle components (PC) (multivariate01, ) of selected variables and then selecting important PC factors from regression with norm.
Let us assume, without loss of generality, that the first dimensions of the vector are nonzero and all the rest are zeros. Selecting the first corresponding entries from the data vector leads to a modified data vector denoted as
(9) 
where denotes the indices where . Next, we perform principle component regression of market arrival series on selected variables with penalty as shown in Eq. 10.
(10) 
where is the matrix obtained from stacking which is the projection of on the principal component. Collecting different for different time instances we create a target vector .
After this exercise, we are left with very few factors . Once we select important PC factors, we run linear regression with arrivals as dependent variable and the PC factors as predictors as shown in Eq.11. Note that ’s need not be positive as production is postively related to NDVI only certain part of lifecycle of crop. e.g. during harvesting period decline in NDVI followed by large arrival in market. Also, factors associated with water body are negative related, as low NDVI in such regions means good water availability, which is good for crop.
(11) 
In all above equations, we took as log arrival due to the constraint that arrival cannot be less than zero and sensitivity at high and low scale of is not same with respect to any changes in factors. We validated these assumptions and ran the experiments for several major Tur markets in Karnataka.
We compare our proposed model, denoted as RegPCR, against standard tree based methods, RandomForest and Boosting, time series based method ARIMA (Hyndman08automatictime, ; timeseries01, )
and other methods such as ridge and principle component regression. Random forest is known to take care of high dimensional data without overfitting and with relatively minimum effort in parameter tunning
(randomforest, ; Statistics01randomforests, ). Gradient Boosting (GB)
(Friedman2002, ; Ridgeway05, ) is known to perform better than other tree based methods, but tends to overfit with increasing number of trees. We experimented with GB’s parameters (no. of trees, shrinkage and tree depth) and selected best combination. In ARIMA, a univariate arrival time series is modeled and used for prediction for next month. We compute the mean absolute error (MAE) between the predicted arrivals and actual values obtained from the market (Eq. 12).(12) 
The comparative results are shown in Fig. 4. Our method performs better than other methods in most markets and is also the clear winner for the aggregated performance. Tree based methods  Random forest, Gradient Boosting do not perform well across most markets. It is due to the fact that with such a small sample size computation intensive methods often fail. We also evaluated other parametric methods appropriate in high dimension data like ours. Ridge regression not consistent across markets but perform better than RandomForest and Boosting. We also observed that direct application of Principle component regression perform poorly across markets. This could be due to noise contribution by nonrelevant locations.
We explore relation among different market arrivals with total arrival of the commodity in the state. The primary hypothesis is the notion that primary markets within a state are strong indicators of the total arrival to the state. Consequently, we use the predicted arrivals for the individual markets at Ballary, Kalaburagi, Koppal, and Raichur to predict the total arrival of Tur for Karnataka.
(13) 
where . The simple regression model for total arrival in Karnataka , gives an adjustedRsquare with a pvalue .
4. Data driven Insights from Market Analysis
4.1. Location Identification
We analyzed the models presented in Sec. 3, to derive useful insights from the commodity arrival prediction workflow. For each market, we selected important variables (latitudelongitude pairs) identified by our model. Since we predict for multiple times, we accumulate the importance of a region based on how important it was for one particular prediction (regression coefficient magnitude) and also weight it with the number of times this location was dominant across the various predictions. We selected important variables which are consistently picked by elasticnet regression and have large absolute regression coefficients (Eq. 3, Eq. 9).
We take two markets to explain our findings. Kalaburgi, the strongest Tur market and Raichur  relatively weak market which tends to grow more than one crop in same season. For Kalaburgi, the important locations are distributed over entire region, thereby pointing to the fact that all the locations are similar in their importance. On the other hand Raichur provide an interesting picture where we find multiple clusters  some are around water bodies and other are large arable land. Note that in Raichur important locations are restricted to only some part of region, highlighting locations which are contributing in Tur production. It also implies other arable locations which are not selected as important contribute in different agricommodity production, not in Tur. The locations with their relative importance (circle size and color for disambiguation) are shown in Fig. 5.
Our method predicts strong region association with the locations which are close to water bodies. This pattern is replicated across the other locations as well, but fades when the market is not among top producer of the crop itself. More explorations along these lines has been marked as a future research activity.
4.2. CropCutting Experiments
One of the important problems which the government in India is facing is to come up with better locations for the CropCutting Experiments (CCEs) (CCE, ). The current guidelines for performing CCEs have the following salient features

Crop Cutting Experiments (CCE) under scientifically designed General Crop Estimation Surveys (GCES).

Around 950 thousand CCEs across India.

Stratified multistage random sampling: Tehsil / Taluk Revenue Village Survey Number / Field Experimental Plot (Specified size / shape)

80120 CCEs for a crop in a major district.
Although the strategy to perform the CCEs looks well motivated from the salient features, the survey results done later tell a different story (KAConf, ). According to the published work, the CCE estimates for Tur, performed by the central government agency and the state government agencies differed from to ! Many of these discrepancies have been attributed to the sampling system which can be also influenced by external factors. Our method can identify the important locations for CCEs for the next cycle based on historical evidence and remote sensing data, thereby ensuring transparency as well as equitable gains. Additionally, identified farm lands in important locations can also be used as investigating points to identify real reason (soil quality, choice of fertilizer, farm management practices) for their consistent performance and replicating it as best practice across regions.
4.3. Price Prediction
Agriculture commodities arrive to the market and are sold to buyers on various days in a month. For each day, the government records the mode, min and max price for each commodity over each market. There are some differences in price distribution from market to market, and the individual daily distributions across the markets are noisy as well. To make the price predictions useful for the state government’s policy making efforts, we propose to predict the monthly average price for the state instead of predicting daily variations for each market. We use the predicted arrivals to model the commodity price in the market.
It is a well known principle of macroeconomics that as the commodity arrival increases / decreases price goes down / up respectively. Generally, a convex function describes the relationship between price (P) and supply/arrivals (A), given by
(14) 
where is a nonnegative constant. In Fig. 6, each observation represents monthly arrival and the corresponding price for Tur, for the state of Karnakata. The simple relationship in Eq. 14 does not hold for extended periods of time. It has been observed under various contexts that due to seasonality and many other factors the curve can shift or transform in various ways.
A closer look at Fig. 6, shows that there is negative correlation between price and arrival, though price conditioned on arrival has high variance. This is due to the fact that major part of commodity which arrives during the harvest season, is consumed all over the year, and therefore, high surge in arrivals may not reflect proportional drop in price and viceversa. The price is also affected by irregular interventions by the government either by declaring changes in minimum support prices to support farmers in case of over supply or by announcing import initiative or export ban to cool off inflation due to crop failure.
Let be price time series and is log arrival time series. To incorporate consumption effect we take exponential weighted average on past observations. We define modified arrival at time t as
(15) 
(16) 
where
(17) 
such that is the exponential decay factor and the index runs over last twelve months. Next, we take differences with duration , and define
(18) 
For , the above equation calculates extra commodity arrived at month t compared to last year. This additional commodity should have effect on the price variations observed for the current and coming months. We model =  as a function of , where d=12. The variation of over the weighted arrivals is shown in Fig. 7.
Once we predict arrival for the next month , we can predict the corresponding price for the next month and for a few more months in the future. Let us assume we want to predict the difference in price, , based on realized arrivals available till time , and the estimated arrival . Eq. 19 denotes this prediction task for the next three months indexed by .
(19) 
We observed that for near term predictions our method works well, possibly because it captures monthly price movement and also takes care of annual commodity availability. The comparison of our model against ARIMA model over last 12 months is shown in Fig. 8. As we move from predicting the prices for the next month, to upto 3 months in the future, the errors by both the models increase.
5. Conclusion
In this paper, we propose a regression framework to predict arrival quantity of Tur to the major markets of Karnataka, India. The proposed method works for extremely unbalanced data setup, where there are less than a hundred observations, but the dimension of those observations can range in millions.
Our entire approach was motivated by the fact that developing nations do not have large amounts of clean recorded data, but these can still be used in principled frameworks to derive value out of them. We use the extremely limited market data and infuse it with large amounts of remote sensing data to arrive at the unbalanced system of equations which is then solved to predict market factors. We utilize the predicted commodity arrivals to draw important insights which can potentially be used by the policy makers to design better optimized cropcutting experiments and modify farm practices by identifying important attributes for consistent performance. We were able to use the predicted arrivals to forecast the corresponding price outlook for three months in the future. Our results show that we outperform traditional methods like Random Forest, Boosting, Ridge Regression, ARIMA and Principle component regression by considerable margins. We ruled out the application of neural techniques due to the extreme scarcity of data. We propose to investigate price prediction in more detail as it is an important problem by itself. We also propose to continue studying the key factors found in final model and explore HSI based analysis.
6. Acknowldegment
We would like to acknowldege Niranjan Nayak, Prashant Gupta, Sushil Chordia and KAPC chairman T.N. Prakash Kammardi for their guidance and support for this project.
References
 (1) D. Donaldson and A. Storeygard, “The view from above: Applications of satellite data in economics,” Journal of Economic Perspectives, vol. 30, no. 4, pp. 171–198, 2016.
 (2) https://modis.gsfc.nasa.gov/.
 (3) https://en.wikipedia.org/wiki/Moderateresolution_imaging_spectroradiometer.
 (4) M. A. D’Iorio and J. Cihlar, “Relationship between AVHRR NDVI and environmental parameters,” in 12th Canadian Symposium on Remote Sensing Geoscience and Remote Sensing Symposium,, vol. 3, July 1989, pp. 1326–1330.
 (5) N. S. Novella and W. M. Thiaw, “African rainfall climatology version 2 for famine early warning systems,” Journal of Applied Meteorology and Climatology, vol. 52, pp. 588–606, 03 2013.

(6)
F. Rembold, C. Atzberger, I. Savin, and O. Rojas, “Using low resolution satellite imagery for yield prediction and yield anomaly detection,”
Remote Sensing, vol. 5, no. 4, pp. 1704–1733, 2013.  (7) P. Doraiswamy, J. Hatfield, T. Jackson, B. Akhmedov, J. Prueger, and A. Stern, “Crop condition and yield simulations using landsat and modis,” Remote Sensing of environment, vol. 92, no. 4, pp. 548–559, 2004.
 (8) B. Kamble, A. Kilic, and K. Hubbard, “Estimating crop coefficients using remote sensingbased vegetation index,” Remote Sensing, vol. 5, no. 4, pp. 1588–1602, 2013.
 (9) T. Aviv and V. LundsgaardNielsen, “Ensemble of cubist models for soy yield prediction using soil features and remote sensing variables,” Data Science for Intelligent Food, Energy and Water, KDD Workshop, 2017.
 (10) J. You, X. Li, M. Low, D. Lobell, and S. Ermon, “Deep gaussian process for crop yield prediction based on remote sensing data,” in AAAI. AAAI Press, 2017, pp. 4559–4566.
 (11) http://krishimaratavahini.kar.nic.in/.
 (12) http://eands.dacnet.nic.in/PDF/MSPkharif201617.pdf.
 (13) H. Zou and T. Hastie, “Regularization and variable selection via the elastic net,” Journal of the Royal Statistical Society, Series B, vol. 67, pp. 301–320, 2005.
 (14) J. Friedman, T. Hastie, and R. Tibshirani, “Regularization paths for generalized linear models via coordinate descent,” Journal of Statistical Software, 2009.
 (15) J. F. Trevor Hastie, Robert Tibshirani, The Elements of Statistical Learning: Data Mining, Inference, and Prediction. Springer, 2003.
 (16) J. Friedman, T. Hastie, H. Höfling, and R. Tibshirani, “Pathwise coordinate optimization,” Annals of Applied Statistics, Tech. Rep., 2007.
 (17) T. Anderson, An Introduction to Multivariate Statistical Analysis, ser. Wiley Series in Probability and Statistics. Wiley, 2003.
 (18) R. J. Hyndman and Y. Kh, “Automatic time series forecasting: The forecast package for r,” Journal of Statistical Software, 2008.
 (19) P. Brockwell and R. Davis, Introduction to Time Series and Forecasting, ser. Introduction to Time Series and Forecasting. Springer, 2002.
 (20) A. Liaw and M. Wiener, “Classification and Regression by randomForest,” R News, vol. 2, no. 3, pp. 18–22, 2002.
 (21) L. Breiman, “Random forests,” Mach. Learn., vol. 45, no. 1, pp. 5–32, Oct. 2001.
 (22) J. H. Friedman, “Stochastic gradient boosting,” Comput. Stat. Data Anal., vol. 38, no. 4, pp. 367–378, Feb. 2002.
 (23) G. Ridgeway, “Generalized boosted models: A guide to the gbm package,” 2005.
 (24) https://ec.europa.eu/agriculture/sites/agriculture/files/events/2015/globcast/10neetu_en.pdf.
 (25) http://des.kar.nic.in/sites/Conference%20Files/2%20Consolidated%20%20All%20Articles%20%20Agriculture1.pdf.