<jupyter_start><jupyter_text>CPSC 330 hw7<jupyter_code>import numpy as np
import pandas as pd
### BEGIN SOLUTION
from sklearn.impute import SimpleImputer
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler, OrdinalEncoder, OneHotEncoder

from sklearn.linear_model import Ridge
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import GridSearchCV, RandomizedSearchCV

from sklearn.metrics import r2_score
### END SOLUTION<jupyter_output><empty_output><jupyter_text>Instructionsrubric={points:5}Follow the [homework submission instructions](https://github.students.cs.ubc.ca/cpsc330-2019w-t2/home/blob/master/docs/homework_instructions.md).   Exercise 1: time series predictionIn this exercise we'll be looking at a [dataset of avocado prices](https://www.kaggle.com/neuromusic/avocado-prices). You should start by downloading the dataset. As usual, please do not commit it to your repos.<jupyter_code>df = pd.read_csv("avocado.csv", parse_dates=["Date"], index_col=0)
df.head()
df.shape
df["Date"].min()
df["Date"].max()<jupyter_output><empty_output><jupyter_text>It looks like the data ranges from the start of 2015 to March 2018 (~2 years ago), for a total of 3.25 years or so. Let's split the data so that we have a 6 months of test data.<jupyter_code>split_date = '20170925'
df_train = df[df["Date"] <= split_date]
df_test  = df[df["Date"] >  split_date]
assert len(df_train) + len(df_test) == len(df)<jupyter_output><empty_output><jupyter_text>1(a)rubric={points:3}In the Rain is Australia dataset from Lecture 16, we had different measurements for each Location. What about this dataset: for which categorical feature(s), if any, do we have separate measurements? Justify your answer by referencing the dataset.  BEGIN SOLUTION<jupyter_code>df.sort_values(by="Date").head()<jupyter_output><empty_output><jupyter_text>From the above, we definitely see measurements on the same day at different regresion. Let's now group by region.<jupyter_code>df.sort_values(by=["region", "Date"]).head()<jupyter_output><empty_output><jupyter_text>From the above we see that, even in Albany, we have two measurements on the same date. This seems to be due to the type of avocado.<jupyter_code>df.sort_values(by=["region", "type", "Date"]).head()<jupyter_output><empty_output><jupyter_text>Great, now we have a sequence of dates with a single row per date. So, the answer is that we have a separate timeseries for each combination of `region` and `type`.  END SOLUTION  1(b)rubric={points:3}In the Rain in Australia dataset, the measurements were generally equally spaced but with some exceptions. How about with this dataset? Justify your answer by referencing the dataset.  BEGIN SOLUTION I think it's not unreasonable to do this on `df` rather than `df_train`, but either way is fine.<jupyter_code>for name, group in df.groupby(['region', 'type']):
    print("%-40s %s" % (name, group["Date"].sort_values().diff().min()))
for name, group in df.groupby(['region', 'type']):
    print("%-40s %s" % (name, group["Date"].sort_values().diff().max()))<jupyter_output>('Albany', 'conventional')               7 days 00:00:00
('Albany', 'organic')                    7 days 00:00:00
('Atlanta', 'conventional')              7 days 00:00:00
('Atlanta', 'organic')                   7 days 00:00:00
('BaltimoreWashington', 'conventional')  7 days 00:00:00
('BaltimoreWashington', 'organic')       7 days 00:00:00
('Boise', 'conventional')                7 days 00:00:00
('Boise', 'organic')                     7 days 00:00:00
('Boston', 'conventional')               7 days 00:00:00
('Boston', 'organic')                    7 days 00:00:00
('BuffaloRochester', 'conventional')     7 days 00:00:00
('BuffaloRochester', 'organic')          7 days 00:00:00
('California', 'conventional')           7 days 00:00:00
('California', 'organic')                7 days 00:00:00
('Charlotte', 'conventional')            7 days 00:00:00
('Charlotte', 'organic')                 7 days 00:00:00
('Chicago', 'conventional')              7 days 00:00:00
('Chicago', 'organic')         [...]<jupyter_text>It looks almost perfect - just organic avocados in WestTexNewMexico seems to be missing a couple measurements.<jupyter_code>name
group["Date"].sort_values().diff().value_counts()<jupyter_output><empty_output><jupyter_text>So, in one case there's a 2-week jump, and in one cast there's a 3-week jump.<jupyter_code>group["Date"].sort_values().reset_index(drop=True).diff().sort_values()<jupyter_output><empty_output><jupyter_text>We can see the anomalies occur at index 48 and 127. (Note: I had to `reset_index` because the index was not unique to each row.)<jupyter_code>group["Date"].sort_values().reset_index(drop=True)[45:50]<jupyter_output><empty_output><jupyter_text>We can spot the first anomaly: a 2-week jump from Nov 29, 2015 to Dec 13, 2015.<jupyter_code>group["Date"].sort_values().reset_index(drop=True)[125:130]<jupyter_output><empty_output><jupyter_text>And we can spot the second anomaly: a 3-week jump from June 11, 2017 to July 2, 2017.  END SOLUTION  1(c)rubric={points:1}In the Rain is Australia dataset, each location was a different place in Australia. For this dataset, look at the names of the regions. Do you think the regions are also all distinct, or are there overlapping regions? Justify your answer by referencing the data.  BEGIN SOLUTION<jupyter_code>df["region"].unique()<jupyter_output><empty_output><jupyter_text>There seems to be a hierarchical structure here: `TotalUS` is split into bigger regions like `West`, `Southeast`, `Northeast`, `Midsouth`; and `California` is split into cities like `Sacramento`, `SanDiego`, `LosAngeles`. It's a bit hard to figure out what's going on.<jupyter_code>df.query("region == 'TotalUS' and type == 'conventional' and Date == '20150104'")["Total Volume"].values[0]
df.query("region != 'TotalUS' and type == 'conventional' and Date == '20150104'")["Total Volume"].sum()<jupyter_output><empty_output><jupyter_text>Since the individual regions sum up to more than the total US, it seems that some of the other regions are double-counted, which is consistent with a hierarchical structure. For example, Los Angeles is probalby double counted because it's within `LosAngeles` but also within `California`. What a mess!  END SOLUTION We will use the entire dataset despite any location-based weirdness uncovered in the previous part.We will be trying to forecast the avocado price, which is the `AveragePrice` column. The function below is adapted from Lecture 16, with some improvements.<jupyter_code>def create_lag_feature(df, orig_feature, lag, groupby, new_feature_name=None, clip=False):
    """
    Creates a new feature that's a lagged version of an existing one.
    
    NOTE: assumes df is already sorted by the time columns and has unique indices.
    
    Parameters
    ----------
    df : pandas.core.frame.DataFrame
        The dataset.
    orig_feature : str
        The column name of the feature we're copying
    lag : int
        The lag; negative lag means values from the past, positive lag means values from the future
    groupby : list
        Column(s) to group by in case df contains multiple time series
    new_feature_name : str
        Override the default name of the newly created column
    clip : bool
        If True, remove rows with a NaN values for the new feature
    
    Returns
    -------
    pandas.core.frame.DataFrame
        A new dataframe with the additional column added.
    """
        
    if new_feature_name is None:
        if lag < 0:
            new_feature_name = "%s_lag%d" % (orig_feature, -lag)
        else:
            new_feature_name = "%s_ahead%d" % (orig_feature, lag)
    
    new_df = df.assign(**{new_feature_name : np.nan})
    for name, group in new_df.groupby(groupby):        
        if lag < 0: # take values from the past
            new_df.loc[group.index[-lag:],new_feature_name] = group.iloc[:lag][orig_feature].values
        else:       # take values from the future
            new_df.loc[group.index[:-lag], new_feature_name] = group.iloc[lag:][orig_feature].values
            
    if clip:
        new_df = new_df.dropna(subset=[new_feature_name])
        
    return new_df<jupyter_output><empty_output><jupyter_text>We first sort our dataframe properly:<jupyter_code>df_sort = df.sort_values(by=["region", "type", "Date"]).reset_index(drop=True)
df_sort<jupyter_output><empty_output><jupyter_text>We then call `create_lag_feature`. This creates a new column in the dataset `AveragePriceNextWeek`, which is the following week's `AveragePrice`. We have set `clip=True` which means it will remove rows where the target would be missing.<jupyter_code>df_hastarget = create_lag_feature(df_sort, "AveragePrice", +1, ["region", "type"], "AveragePriceNextWeek", clip=True)
df_hastarget<jupyter_output><empty_output><jupyter_text>I will now split the data:<jupyter_code>df_train = df_hastarget[df_hastarget["Date"] <= split_date]
df_test  = df_hastarget[df_hastarget["Date"] >  split_date]<jupyter_output><empty_output><jupyter_text>1(d)rubric={points:1}Why was it reasonable for me to do this operation _before_ splitting the data, despite the fact that this usually constitutes a violation of the Golden Rule?  BEGIN SOLUTIONBecause we were only looking at the dates and creating the future feature. The difference is that the very last time point in our training set now contains the average price from the first time point in our test set. This is a realistic scenario if we wre actually using this model to forecast, so it's not a major concern. END SOLUTION  1(e)rubric={points:1}Next we will want to build some models to forecast the average avocado price a week in advance. Before we start with any ML, let's try a baseline: just predicting the previous week's `AveragePrice`. What $R^2$ do you get with this approach?  BEGIN SOLUTION<jupyter_code>r2_score(df_train["AveragePriceNextWeek"], df_train["AveragePrice"])
r2_score(df_test["AveragePriceNextWeek"], df_test["AveragePrice"])<jupyter_output><empty_output><jupyter_text>Interesting that this is a less effective prediction strategy in the later part of the dataset. I guess that means the price was fluctuating more in late 2017 / early 2018?  END SOLUTION  1(f)rubric={points:10}Build some models to forecast the average avocado price. Experiment with a few approachs for encoding the date. Justify the decisions you make. Which approach worked best? Report your test score and briefly discuss your results.Benchmark: you should be able to achieve $R^2$ of at least 0.79 on the test set. I got to 0.80, but not beyond that. Let me know if you do better!Note: because we only have 2 splits here, we need to be a bit wary of overfitting on the test set. Try not to test on it a ridiculous number of times. If you are interested in some proper ways of dealing with this, see for example sklearn's [TimeSeriesSplit](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.TimeSeriesSplit.html), which is like cross-validation for time series data.  BEGIN SOLUTION<jupyter_code>df_train.head()
(df_train.loc[:, "Small Bags": "XLarge Bags"].sum(axis=1) - df_train["Total Bags"]).abs().max()<jupyter_output><empty_output><jupyter_text>It seems that `Total Bags` is (approximately) the sum of the other 3 bag features, so I will drop `Total Bags`.<jupyter_code>(df_train.loc[:, "4046": "4770"].sum(axis=1) - df_train["Total Volume"]).abs().max()<jupyter_output><empty_output><jupyter_text>It seems that `Total Volume` is _not_ the sum of the 3 avocado types, so I will keep all 4 columns.<jupyter_code>df_train.info()<jupyter_output><class 'pandas.core.frame.DataFrame'>
Int64Index: 15441 entries, 0 to 18222
Data columns (total 14 columns):
 #   Column                Non-Null Count  Dtype         
---  ------                --------------  -----         
 0   Date                  15441 non-null  datetime64[ns]
 1   AveragePrice          15441 non-null  float64       
 2   Total Volume          15441 non-null  float64       
 3   4046                  15441 non-null  float64       
 4   4225                  15441 non-null  float64       
 5   4770                  15441 non-null  float64       
 6   Total Bags            15441 non-null  float64       
 7   Small Bags            15441 non-null  float64       
 8   Large Bags            15441 non-null  float64       
 9   XLarge Bags           15441 non-null  float64       
 10  type                  15441 non-null  object        
 11  year                  15441 non-null  int64         
 12  region                15441 non-null  object        
 13  AveragePriceNext[...]<jupyter_text>It seems there are no null values, so I will not do any imputation. Will plot a single time series for exploration purposes:<jupyter_code>df_train.query("region == 'TotalUS'").set_index("Date").groupby("type")["AveragePrice"].plot(legend=True);
df_train.query("region == 'TotalUS' and type == 'conventional'").plot(x="Date", y="Total Volume");<jupyter_output><empty_output><jupyter_text>We see some seasonality in the total volume, but not much in the average price - interesting. I will not scale the `AveragePrice` because I am not scaling `AveragePriceNextWeek` either, and it may be helpful to keep them the same. Alternatively, it may have been effective to predict the _change_ in price instead of next's week's price.<jupyter_code>numeric_features     = ["Total Volume", "4046", "4225", "4770", "Small Bags", "Large Bags", "XLarge Bags", "year"]
categorical_features = ["type", "region"]
keep_features        = ["AveragePrice"]
drop_features        = ["Date", "Total Bags"] 
target_feature       = "AveragePriceNextWeek"<jupyter_output><empty_output><jupyter_text>Next, I grab the `preprocess_features` function from Lecture 16, with a minor modification to allow un-transformed features via `keep_features`:<jupyter_code>def preprocess_features(df_train, df_test, 
                        numeric_features, 
                        categorical_features, 
                        keep_features, 
                        drop_features, 
                        target_feature):
    
    all_features = numeric_features + categorical_features + keep_features + drop_features + [target_feature]
    if set(df_train.columns) != set(all_features):
        print("Missing columns", set(df_train.columns) - set(all_features))
        print("Extra columns", set(all_features) - set(df_train.columns))
        raise Exception("Columns do not match")
        
    # Put the columns in the order we want
    df_train = df_train[all_features]
    df_test = df_test[all_features]
    
    numeric_transformer = Pipeline([
        ('imputer', SimpleImputer(strategy='median')),
        ('scaler', StandardScaler())
    ])
    categorical_transformer = Pipeline([
        ('imputer', SimpleImputer(strategy='most_frequent')),
        ('onehot', OneHotEncoder(sparse=False, drop='first'))
    ])
    preprocessor = ColumnTransformer([
        ('numeric', numeric_transformer, numeric_features),
        ('categorical', categorical_transformer, categorical_features)
    ], remainder='passthrough')
    preprocessor.fit(df_train);

    if len(categorical_features) > 0:
        ohe = preprocessor.named_transformers_['categorical'].named_steps['onehot']
        ohe_feature_names = list(ohe.get_feature_names(categorical_features))
        new_columns = numeric_features + ohe_feature_names + keep_features + drop_features + [target_feature]
    else:
        new_columns = all_features

    X_train_enc = pd.DataFrame(preprocessor.transform(df_train), index=df_train.index, columns=new_columns)
    X_test_enc  = pd.DataFrame(preprocessor.transform(df_test),  index=df_test.index,  columns=new_columns)
    
    X_train_enc = X_train_enc.drop(columns=drop_features + [target_feature])
    X_test_enc  = X_test_enc.drop( columns=drop_features + [target_feature])
    
    y_train = df_train[target_feature]
    y_test  = df_test[ target_feature]
    
    return X_train_enc, y_train, X_test_enc, y_test
df_train_enc, y_train, df_test_enc, y_test = preprocess_features(df_train, df_test, 
                        numeric_features, 
                        categorical_features, 
                        keep_features, 
                        drop_features, 
                        target_feature)
df_train_enc.head()
lr = Ridge()
lr.fit(df_train_enc, y_train);
lr.score(df_train_enc, y_train)
lr.score(df_test_enc, y_test)
lr_coef = pd.DataFrame(data=np.squeeze(lr.coef_), index=df_train_enc.columns, columns=["Coef"])
lr_coef.sort_values(by="Coef", ascending=False)<jupyter_output><empty_output><jupyter_text>This is not a very impressive showing. We're doing almost the same as the baseline. Let's see if encoding the date helps at all. We'll try to OHE the month.<jupyter_code>df_train_month = df_train.assign(Month=df_train["Date"].apply(lambda x: x.month))
df_test_month  = df_test.assign( Month=df_test[ "Date"].apply(lambda x: x.month))
df_train_month_enc, y_train, df_test_month_enc, y_test = preprocess_features(df_train_month, df_test_month, 
                        numeric_features, 
                        categorical_features + ["Month"], 
                        keep_features, 
                        drop_features, 
                        target_feature)
df_train_month_enc.head()
lr = Ridge()
lr.fit(df_train_month_enc, y_train);
lr.score(df_train_month_enc, y_train)
lr.score(df_test_month_enc, y_test)<jupyter_output><empty_output><jupyter_text>A tiny bit better.<jupyter_code>pd.DataFrame(data=np.squeeze(lr.coef_), index=df_train_month_enc.columns, columns=["Coef"]).sort_values(by="Coef", ascending=False)<jupyter_output><empty_output><jupyter_text>Let's add some lag features. I'm arbitrarily deciding on 4 lags for `AveragePrice` (the most important feature).<jupyter_code>def add_lags(df):
    df = create_lag_feature(df, "AveragePrice", -1, ["region", "type"])
    df = create_lag_feature(df, "AveragePrice", -2, ["region", "type"])
    df = create_lag_feature(df, "AveragePrice", -3, ["region", "type"])
    df = create_lag_feature(df, "AveragePrice", -4, ["region", "type"])
    return df

df_train_month_lag = add_lags(df_train_month)
df_test_month_lag  = add_lags(df_test_month)
df_train_month_lag
df_train_month_lag_enc, y_train, df_test_month_lag_enc, y_test = preprocess_features(df_train_month_lag, df_test_month_lag, 
                        numeric_features + ["AveragePrice_lag1", "AveragePrice_lag2", "AveragePrice_lag3", "AveragePrice_lag4"], 
                        categorical_features + ["Month"], 
                        keep_features, 
                        drop_features, 
                        target_feature)
lr = Ridge()
lr.fit(df_train_month_lag_enc, y_train);
lr.score(df_train_month_lag_enc, y_train)
lr.score(df_test_month_lag_enc, y_test)<jupyter_output><empty_output><jupyter_text>This did not seem to help.<jupyter_code>pd.DataFrame(data=np.squeeze(lr.coef_), index=df_train_month_lag_enc.columns, columns=["Coef"]).sort_values(by="Coef", ascending=False)<jupyter_output><empty_output><jupyter_text>We can also try a random forest:<jupyter_code>rf = RandomForestRegressor() 
rf.fit(df_train_month_lag_enc, y_train);
rf.score(df_train_month_lag_enc, y_train)
rf.score(df_test_month_lag_enc, y_test)<jupyter_output><empty_output><jupyter_text>For the random forest it may be helpful to model the difference between today and tomorrow. The linear model does not care about this because it just corresponds to changing the coefficient corresponding to `AveragePrice` by 1, but for the random forest it may help:<jupyter_code>rf = RandomForestRegressor() 
rf.fit(df_train_month_lag_enc, y_train - df_train_month_lag_enc["AveragePrice"]);
r2_score(y_train, rf.predict(df_train_month_lag_enc) + df_train_month_lag_enc["AveragePrice"])
r2_score(y_test, rf.predict(df_test_month_lag_enc) + df_test_month_lag_enc["AveragePrice"])<jupyter_output><empty_output><jupyter_text>This massively overfits when we do this shifting. Let's try a simpler model...<jupyter_code>rf = RandomForestRegressor(max_depth=8) 
rf.fit(df_train_month_lag_enc, y_train - df_train_month_lag_enc["AveragePrice"]);
r2_score(y_train, rf.predict(df_train_month_lag_enc) + df_train_month_lag_enc["AveragePrice"])
r2_score(y_test, rf.predict(df_test_month_lag_enc) + df_test_month_lag_enc["AveragePrice"])<jupyter_output><empty_output><jupyter_text>Doesn't realy help. Also, we can just confirm that this shifting has no effect on the linear model (well, a small effect because it's `Ridge` instead of `LinearRegression`, but small):<jupyter_code>lr = Ridge() 
lr.fit(df_train_month_lag_enc, y_train - df_train_month_lag_enc["AveragePrice"]);
r2_score(y_train, lr.predict(df_train_month_lag_enc) + df_train_month_lag_enc["AveragePrice"])
r2_score(y_test, lr.predict(df_test_month_lag_enc) + df_test_month_lag_enc["AveragePrice"])<jupyter_output><empty_output><jupyter_text>Indeed, this is essentially the same score we had before. Overall, adding the month helped, but adding the lagged price was surprisingly unhelpful. Perhaps lagged version of other features would have been better, or other representations of the time of year, or dealing with the regions and avocado types a bit more carefully.   END SOLUTION  1(g)rubric={points:3}We talked a little bit about _seasonality_, which is the idea of a periodic component to the time series. For example, in Lecture 16 we attempted to capture this by encoding the month. Something we didn't discuss is _trends_, which are long-term variations in the quantity of interest. Aside from the effects of climate change, the amount of rain in Australia is likely to vary during the year but less likely to have long-term trends over the years. Avocado prices, on the other hand, could easily exhibit trends: for example avocados may just cost more in 2020 than they did in 2015.Briefly discuss in ~1 paragraph: to what extent, if any, was your model above able to account for seasonality? What about trends?  BEGIN SOLUTIONI tried to take seasonality into account by having the month as an OHE variable. As far as trends are concerned, the year is also a numeric variable in the model, so it could learn that the price in 2017 is higher than in 2015, say. However, there are very few years in the training set (2015, 16, 17), so that is not a lot of data to learn from. Perhaps including the number of months since the start of the dataset, or something like that, would enable the model to do a bit better with trends. Nonetheless, extrapolating is very hard so we can't necessarily trust our models' handing of trend.<jupyter_code>pd.DataFrame(data=np.squeeze(lr.coef_), index=df_train_month_lag_enc.columns, columns=["Coef"]).loc["year"]<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09778904133531417, 'fraction_numerical': 0.020274403565498558, 'mean_word_length': 4.244500458295142, 'pattern_counts': {'":': 2, '<': 136, '<?xml version=': 0, '>': 135, 'https://': 3, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '33', 'dir': 'jupyter-structured-clean-dedup', 'id': '8570', 'license': 'MIT', 'n_tokens_mistral': 7435, 'n_tokens_neox': 6859, 'n_words': 2169, 'path': 'junk/.ipynb_checkpoints/hw7-checkpoint.ipynb', 'repo_name': 'Lavendulaa/programming-in-python-for-data-science', 'size': '22886'}