{"cells":[{"cell_type":"markdown","id":"2e980a87-afac-46e1-b8e5-758c16be989d","metadata":{},"outputs":[],"source":["

\n"," \n"," \"Skills\n"," \n","

\n"]},{"cell_type":"markdown","id":"d6d03fff-24fa-4a06-a4d7-307dc0272c94","metadata":{},"outputs":[],"source":["## NOAA Weather Data Analysis - JFK Airport (New York)\n","\n","This notebook focuses on analyzing and forecasting weather patterns using the **NOAA Weather Dataset** collected from **JFK Airport in New York**. The dataset comprises **114,546 hourly observations** of **12 key climatological variables**, including temperature, wind speed, humidity, and pressure.\n","This notebook teaches the user to extract, clean and analyze sample weather data and predict weather trends to help airports schedule better flight times. \n","\n","\n","The notebook is organized into three main parts:\n","\n","---\n","\n","## Part 1: Data Cleaning\n","\n","In this section, we prepare the raw data for analysis by:\n","\n","- Removing unnecessary or redundant columns to retain only relevant numerical features\n","- Converting data types and cleaning inconsistencies\n","- Handling missing values with appropriate filling strategies\n","- Encoding categorical weather features for downstream analysis\n","\n","---\n","\n","## Part 2: Exploratory Data Analysis (EDA)\n","\n","Here, we perform visual and statistical exploration of the cleaned dataset:\n","\n","- Load the cleaned data\n","- Generate insightful visualizations of key variables\n","- Identify trends, patterns, and seasonal effects in the time-series data\n","\n","---\n","\n","## Part 3: Time Series Forecasting\n","\n","This section focuses on predicting future temperatures using time-series models:\n","\n","- Load the cleaned and preprocessed data\n","- Establish baseline forecasting models\n","- Train and evaluate advanced statistical forecasting techniques\n"]},{"cell_type":"markdown","id":"cf471622-4eb2-4665-9c6e-224bb8960867","metadata":{},"outputs":[],"source":["## Part 1: Data Cleaning\n","\n"]},{"cell_type":"code","id":"48cf287a-b77d-4bcc-af7f-f1db394772d4","metadata":{},"outputs":[],"source":["%%capture\n!pip install pandas\n#Import and configure the required modules.\n# Define required imports\nimport pandas as pd\nimport numpy as np\nimport sys\nimport re\n# These set pandas max column and row display in the notebook\npd.set_option('display.max_columns', 50)\npd.set_option('display.max_rows', 50)"]},{"cell_type":"markdown","id":"39a2b611-4723-47d6-a625-39d7db5a6ec1","metadata":{},"outputs":[],"source":["### Read the Raw Data\n","We start by reading in the raw dataset, displaying the first few rows of the dataframe, and taking a look at the columns and column types present.\n"]},{"cell_type":"code","id":"04b6e057-779b-4bd7-84ad-b47c5cb39c96","metadata":{},"outputs":[],"source":["%%capture\n# Using pandas to read the data \n# Since the `DATE` column consists date-time information, we use Pandas parse_dates keyword for easier data processing\n\nraw_data = pd.read_csv(\"https://cf-courses-data.s3.us.cloud-object-storage.appdomain.cloud/fJKLoMKtgpDAl5MkA6rf2A/jfk-weather.csv\", parse_dates=['DATE'])\nraw_data.head()"]},{"cell_type":"code","id":"d5d89ad7-a5ed-405b-a1a7-29def5c81523","metadata":{},"outputs":[],"source":["raw_data.dtypes"]},{"cell_type":"markdown","id":"b418ace8-e763-47a1-b1b5-6063c4711789","metadata":{},"outputs":[],"source":["## Clean the Data\n","As you can see above, there are a lot of fields which are non-numerical - usually these will be fields that contain text or categorical data, e.g. HOURLYSKYCONDITIONS.\n","\n","There are also fields - such as the main temperature field of interest HOURLYDRYBULBTEMPF - that we expect to be numerical, but are instead object type. This often indicates that there may be missing (or null) values, or some other unusual readings that we may have to deal with (since otherwise the field would have been fully parsed as a numerical data type).\n","\n","In addition, some fields relate to hourly observations, while others relate to daily or monthly intervals. For purposes of later exploratory data analysis, we will restrict the dataset to a certain subset of numerical fields that relate to hourly observations.\n","\n","In this section, we refer to the ***NOAA Local Climatological Data Documentation*** to describe the fields and meaning of various values.\n","\n","\n","### Select data columns\n","First, we select only the subset of data columns of interest and inspect the column types.\n"]},{"cell_type":"code","id":"52fec1a1-ec5c-4677-b6a3-0d93f2b3256b","metadata":{},"outputs":[],"source":["# Choose what columns to import from raw data\ncolumn_subset = [\n 'DATE',\n 'HOURLYVISIBILITY',\n 'HOURLYDRYBULBTEMPF',\n 'HOURLYWETBULBTEMPF',\n 'HOURLYDewPointTempF',\n 'HOURLYRelativeHumidity',\n 'HOURLYWindSpeed',\n 'HOURLYWindDirection',\n 'HOURLYStationPressure',\n 'HOURLYPressureTendency',\n 'HOURLYSeaLevelPressure',\n 'HOURLYPrecip',\n 'HOURLYAltimeterSetting'\n]\n\n# Filter dataset to relevant columns\nhourly_data = raw_data[column_subset]\n# Set date index\nhourly_data = hourly_data.set_index(pd.DatetimeIndex(hourly_data['DATE']))\nhourly_data.drop(['DATE'], axis=1, inplace=True)\nhourly_data.replace(to_replace='*', value=np.nan, inplace=True)\nhourly_data.head()"]},{"cell_type":"code","id":"d4e0fcde-d959-4e63-93fd-bfae90e22203","metadata":{},"outputs":[],"source":["#hourly_data.reset_index(inplace=True)"]},{"cell_type":"code","id":"74aa307c-a4c8-415a-9b1a-2ace07e0e37d","metadata":{},"outputs":[],"source":["hourly_data.dtypes"]},{"cell_type":"markdown","id":"fad200ca-2031-4543-bf5b-737bf4ce93ca","metadata":{},"outputs":[],"source":["### Clean up precipitation column\n","From the dataframe preview above, we can see that the column ***HOURLYPrecip*** - which is the hourly measure of precipitation levels - ***contains both NaN and T values***. T specifies trace amounts of precipitation, while NaN means not a number, and is used to denote missing values.\n","\n","We can also inspect the unique values present for the field.\n"]},{"cell_type":"code","id":"ca75dcd9-5bd4-4dcd-8ed1-7d94acad8d9e","metadata":{},"outputs":[],"source":["hourly_data['HOURLYPrecip'].unique()"]},{"cell_type":"markdown","id":"a3713bf8-b981-40d7-84d0-643d57bee0f5","metadata":{},"outputs":[],"source":["We can see that some values end with an s (indicating snow), while there is a strange value 0.020.01s which appears to be an error of some sort. To deal with T values, we will set the observation to be 0. We will also replace the erroneous value 0.020.01s with NaN.\n"]},{"cell_type":"code","id":"b5e96075-9990-4310-a1db-ae60b2573199","metadata":{},"outputs":[],"source":["%%capture\n# Fix imported data\nhourly_data['HOURLYPrecip'].replace(to_replace='T', value='0.00', inplace=True)\nhourly_data['HOURLYPrecip'].replace('0.020.01s', np.nan, inplace=True)"]},{"cell_type":"markdown","id":"55ee8959-2d88-4061-961e-54ee8bf033c5","metadata":{},"outputs":[],"source":["### Convert columns to numerical types\n","Next, we will convert string columns that refer to numerical values to numerical types. For columns such as HOURLYPrecip, we first also drop the non-numerical parts of the value (the s character).\n"]},{"cell_type":"code","id":"7733d72d-eb42-45aa-8e4f-69c9526c5ede","metadata":{},"outputs":[],"source":["# Set of columns to convert\nmessy_columns = column_subset[1:]\n\n# Convert columns to float32 datatype\nfor i in messy_columns:\n hourly_data[i] = hourly_data[i].apply(lambda x: re.sub('[^0-9,.-]', '', x) if type(x) == str else x).replace('', np.nan).astype(('float32'))"]},{"cell_type":"markdown","id":"afc696db-e28a-4948-853f-83c9ae9711cc","metadata":{},"outputs":[],"source":["We can now see that all fields have numerical data type.\n"]},{"cell_type":"code","id":"e28d73f5-cfc8-4fd0-8a0b-9bed127cc983","metadata":{},"outputs":[],"source":["print(hourly_data.info())\nhourly_data.head()"]},{"cell_type":"markdown","id":"0f6377c0-c2ec-4369-a478-dc341131e7d6","metadata":{},"outputs":[],"source":["### Reformat and process data¶\n","Next, we will clean up some of the data columns to ensure their values fall within the parameters defined by the NOAA documentation (referred to above).\n"]},{"cell_type":"code","id":"c6937f91-ea19-45d0-ab0c-c8116cbbdb01","metadata":{},"outputs":[],"source":["# Generate the summary statistics for each column\nhourly_data.describe()"]},{"cell_type":"markdown","id":"626eb79a-d39b-472e-af4d-58b6e1be3387","metadata":{},"outputs":[],"source":["According to the documentation, the HOURLYPressureTendency field should be an integer value in the range [0, 8]. Let's check if this condition holds for this dataset.\n"]},{"cell_type":"code","id":"ef86e2bc-77b4-46cd-b6a9-4eee1056a742","metadata":{},"outputs":[],"source":["# Check if categorical variable HOURLYPressureTendency ever has a non-integer entry outside the bounds of 0-8\ncond = len(hourly_data[~hourly_data['HOURLYPressureTendency'].isin(list(range(0,9)) + [np.nan])])\nprint('Hourly Pressure Tendency should be between 0 and 8: {}'.format(cond == 0))"]},{"cell_type":"markdown","id":"2b5cc1b0-8d4e-463e-a0bc-249c3fa1b0c1","metadata":{},"outputs":[],"source":["Hourly Pressure Tendency should be between 0 and 8: True\n","The HOURLYVISIBILITY should be an integer in the range [0, 10]. Let's check this condition too.\n"]},{"cell_type":"code","id":"99e139ba-a1b7-4e7f-a4fc-f4e8fd326f74","metadata":{},"outputs":[],"source":["# Hourly Visibility should be between 0 and 10\nhourly_data[(hourly_data['HOURLYVISIBILITY'] < 0) | (hourly_data['HOURLYVISIBILITY'] > 10)]"]},{"cell_type":"markdown","id":"aafacf6b-fafc-4351-8628-e792cdf1eb62","metadata":{},"outputs":[],"source":["We find that a couple of observations fall outside the range. These must be spurious data observations and we handle them by replacing them with NaN.\n"]},{"cell_type":"code","id":"5a401dcd-c31c-46ba-a264-1685704d44a0","metadata":{},"outputs":[],"source":["# Replace any hourly visibility figure outside these bounds with nan\nhourly_data.loc[hourly_data['HOURLYVISIBILITY'] > 10, 'HOURLYVISIBILITY'] = np.nan\n\n# Hourly Visibility should be between 0 and 10\ncond = len(hourly_data[(hourly_data['HOURLYVISIBILITY'] < 0) | (hourly_data['HOURLYVISIBILITY'] > 10)])\nprint('Hourly Visibility should be between 0 and 10: {}'.format(cond == 0))"]},{"cell_type":"markdown","id":"8cd8839c-e1f8-41f6-859d-8e8af0ca2435","metadata":{},"outputs":[],"source":["Finally, we check if there are any duplicates with respect to our DATE index and check furthermore that our dates are in the correct order (that is, strictly increasing).\n"]},{"cell_type":"code","id":"ffe1b3b9-504d-490d-bc90-6b87abe4d573","metadata":{},"outputs":[],"source":["cond = len(hourly_data[hourly_data.index.duplicated()].sort_index())\nprint('Date index contains no duplicate entries: {}'.format(cond == 0))"]},{"cell_type":"code","id":"a7e80188-91d0-412e-86cf-755749f386ec","metadata":{},"outputs":[],"source":["# Make sure time index is sorted and increasing\nprint('Date index is strictly increasing: {}'.format(hourly_data.index.is_monotonic_increasing))"]},{"cell_type":"code","id":"c96c7e63-5533-4911-acd1-cffcb224ae69","metadata":{},"outputs":[],"source":["hourly_data.reset_index(inplace=True)"]},{"cell_type":"code","id":"2653eeff-7aa2-4e51-b459-a942cd67f825","metadata":{},"outputs":[],"source":["hourly_data.columns"]},{"cell_type":"markdown","id":"f100d8fe-70a7-482c-b1f7-ec4dc38cc678","metadata":{},"outputs":[],"source":["We will now also replace missing values. For numerical values, we will linearly interpolate between the previous and next valid obvservations. For the categorical HOURLYPressureTendency field, we will replace missing values with the last valid observation.\n"]},{"cell_type":"code","id":"d7406f73-c057-44ce-ba9c-f282d27e3610","metadata":{},"outputs":[],"source":["%%capture\nhourly_data['HOURLYPressureTendency'] = hourly_data['HOURLYPressureTendency'].fillna(method='ffill') # fill with last valid observation\nhourly_data = hourly_data.interpolate(method='linear') # interpolate missing values\nhourly_data.drop(hourly_data.index[0], inplace=True) # drop first row"]},{"cell_type":"code","id":"5ee31b73-d396-485a-b08a-998fa3546425","metadata":{},"outputs":[],"source":["print(hourly_data.info())\nprint()\nhourly_data.head()"]},{"cell_type":"markdown","id":"90f61f50-b246-4b28-a5f6-2015f231204a","metadata":{},"outputs":[],"source":["### Feature encoding\n","The final pre-processing step we will perform will be to handle two of our columns in a special way in order to correctly encode these features. They are:\n","\n","- **HOURLYWindDirection** - wind direction\n","- **HOURLYPressureTendency** - an indicator of pressure changes\n","\n","For HOURLYWindDirection, we encode the raw feature value as two new values, which measure the cyclical nature of wind direction - that is, we are encoding the compass-point nature of wind direction measurements.\n"]},{"cell_type":"code","id":"307cd2fd-cf07-443e-bae1-81d546df6b96","metadata":{},"outputs":[],"source":["# Transform HOURLYWindDirection into a cyclical variable using sin and cos transforms\nhourly_data['HOURLYWindDirectionSin'] = np.sin(hourly_data['HOURLYWindDirection']*(2.*np.pi/360))\nhourly_data['HOURLYWindDirectionCos'] = np.cos(hourly_data['HOURLYWindDirection']*(2.*np.pi/360))\nhourly_data.drop(['HOURLYWindDirection'], axis=1, inplace=True)"]},{"cell_type":"markdown","id":"7faeac47-309c-42cf-8ca7-3ef7a8cf11ab","metadata":{},"outputs":[],"source":["For HOURLYPressureTendency, the feature value is in fact a categorical feature with three levels:\n","\n","- 0-3 indicates an increase in pressure over the previous 3 hours\n","- 4 indicates no change during the previous 3 hours\n","- 5-8 indicates a decrease over the previous 3 hours\n","Hence, we encode this feature into 3 dummy values representing these 3 potential states.\n"]},{"cell_type":"code","id":"d58eda7f-fe8a-49ac-9326-25fcf607f129","metadata":{},"outputs":[],"source":["# Transform HOURLYPressureTendency into 3 dummy variables based on NOAA documentation\nhourly_data['HOURLYPressureTendencyIncr'] = [1.0 if x in [0,1,2,3] else 0.0 for x in hourly_data['HOURLYPressureTendency']] # 0 through 3 indicates an increase in pressure over previous 3 hours\nhourly_data['HOURLYPressureTendencyDecr'] = [1.0 if x in [5,6,7,8] else 0.0 for x in hourly_data['HOURLYPressureTendency']] # 5 through 8 indicates a decrease over previous 3 hours\nhourly_data['HOURLYPressureTendencyConst'] = [1.0 if x == 4 else 0.0 for x in hourly_data['HOURLYPressureTendency']] # 4 indicates no change during previous 3 hours\nhourly_data.drop(['HOURLYPressureTendency'], axis=1, inplace=True)\nhourly_data['HOURLYPressureTendencyIncr'] = hourly_data['HOURLYPressureTendencyIncr'].astype(('float32'))\nhourly_data['HOURLYPressureTendencyDecr'] = hourly_data['HOURLYPressureTendencyDecr'].astype(('float32'))\nhourly_data['HOURLYPressureTendencyConst'] = hourly_data['HOURLYPressureTendencyConst'].astype(('float32'))"]},{"cell_type":"markdown","id":"a1af4c67-ab56-4519-9d8e-72225e37a748","metadata":{},"outputs":[],"source":["### Rename columns\n","Before saving the dataset, we will rename the columns for readability.\n"]},{"cell_type":"code","id":"c8f2a988-7b12-4775-9649-05d9b6e501d0","metadata":{},"outputs":[],"source":["hourly_data.columns"]},{"cell_type":"code","id":"7992b56b-21d5-405e-804b-ed356790c882","metadata":{},"outputs":[],"source":["# define the new column names\ncolumns_new_name = [\n 'DATE',\n 'visibility',\n 'dry_bulb_temp_f',\n 'wet_bulb_temp_f',\n 'dew_point_temp_f',\n 'relative_humidity',\n 'wind_speed',\n 'station_pressure',\n 'sea_level_pressure',\n 'precip',\n 'altimeter_setting',\n 'wind_direction_sin',\n 'wind_direction_cos',\n 'pressure_tendency_incr',\n 'pressure_tendency_decr',\n 'pressure_tendency_const'\n]\n\ncolumns_name_map = {c:columns_new_name[i] for i, c in enumerate(hourly_data.columns)}\n\nhourly_data_renamed = hourly_data.rename(columns=columns_name_map)"]},{"cell_type":"code","id":"1043e2a5-ea53-40c7-b4a5-baae8718989c","metadata":{},"outputs":[],"source":["print(hourly_data_renamed.info())\nprint()\nhourly_data_renamed.head()"]},{"cell_type":"code","id":"3f5e5ee6-183a-446e-876b-4991beee5a1d","metadata":{},"outputs":[],"source":["hourly_data_renamed['DATE'] = pd.to_datetime(hourly_data_renamed['DATE'])\nhourly_data_renamed.set_index('DATE', inplace=True)"]},{"cell_type":"code","id":"4977d02e-7d46-4340-ac84-910bd1380375","metadata":{},"outputs":[],"source":["# Explore some general information about the dataset\nprint('# of megabytes held by dataframe: ' + str(round(sys.getsizeof(hourly_data_renamed) / 1000000,2)))\nprint('# of features: ' + str(hourly_data_renamed.shape[1])) \nprint('# of observations: ' + str(hourly_data_renamed.shape[0]))\nprint('Start date: ' + str(hourly_data_renamed.index[0]))\nprint('End date: ' + str(hourly_data_renamed.index[-1]))\nprint('# of days: ' + str((hourly_data_renamed.index[-1] - hourly_data_renamed.index[0]).days))\nprint('# of months: ' + str(round((hourly_data_renamed.index[-1] - hourly_data_renamed.index[0]).days/30,2)))\nprint('# of years: ' + str(round((hourly_data_renamed.index[-1] - hourly_data_renamed.index[0]).days/365,2)))"]},{"cell_type":"markdown","id":"8b99da06-a64a-48e7-8956-e9e44eec984a","metadata":{},"outputs":[],"source":["## Save the Cleaned Data\n","Finally, we save the cleaned dataset.\n"]},{"cell_type":"code","id":"46ce2239-fd42-479b-86cf-55aa497a8b1c","metadata":{},"outputs":[],"source":["hourly_data_renamed.to_csv(\"jfk_weather_cleaned.csv\", float_format='%g')"]},{"cell_type":"code","id":"4a981709-8120-4d4c-894f-82dda8cf0d34","metadata":{},"outputs":[],"source":["hourly_data_renamed"]},{"cell_type":"markdown","id":"634b0e38-c1a8-46e6-a59b-21b1f1cd5c0d","metadata":{},"outputs":[],"source":["## Part 2: Exploratory Data Analysis (EDA)\n"]},{"cell_type":"code","id":"732db809-a452-43f5-b13a-4d37ed0e4a34","metadata":{},"outputs":[],"source":["%%capture\n# Installing packages needed for data processing and visualization\n!pip install pandas matplotlib seaborn numpy "]},{"cell_type":"code","id":"7a9e6fe8-2554-4c42-aa4f-44cdabc7012d","metadata":{},"outputs":[],"source":["import seaborn as sns\nfrom pandas import DataFrame as df\nfrom matplotlib import pyplot as plt\nplt.rcParams['figure.dpi'] = 160"]},{"cell_type":"code","id":"4be40e3f-241f-4ee3-8250-2c59b2b25c11","metadata":{},"outputs":[],"source":["data = pd.read_csv(\"jfk_weather_cleaned.csv\")"]},{"cell_type":"markdown","id":"5b31a50f-ce45-4447-8f97-7d7698746614","metadata":{},"outputs":[],"source":["### Read the Cleaned Data¶\n","We start by reading in the cleaned dataset that was saved in ***Part 1 - Data Cleaning.**\n"]},{"cell_type":"code","id":"40c75346-57a2-4fe4-a65b-3a9ebf6d1b55","metadata":{},"outputs":[],"source":["data.head()"]},{"cell_type":"code","id":"660667e7-50bc-48b4-a7b4-97d90af12f80","metadata":{},"outputs":[],"source":["data = pd.read_csv(\"jfk_weather_cleaned.csv\", parse_dates=['DATE'])\n# Set date index\ndata = data.set_index(pd.DatetimeIndex(data['DATE']))\ndata.drop(['DATE'], axis=1, inplace=True)\ndata.head()"]},{"cell_type":"markdown","id":"5dc953f7-0b1b-44f1-a2f7-2c93be7a70e2","metadata":{},"outputs":[],"source":["## Visualize the Data\n","In this section we visualize a few sections of the data, using matplotlib's pyplot module.\n"]},{"cell_type":"code","id":"8bb88629-3dd8-4c0d-9d78-044a6cd7435e","metadata":{},"outputs":[],"source":["# Columns to visualize\nplot_cols = ['dry_bulb_temp_f', 'relative_humidity', 'wind_speed', 'station_pressure', 'precip']"]},{"cell_type":"markdown","id":"78d66ee3-66c8-4204-bc8e-fb6305591a26","metadata":{},"outputs":[],"source":["#### Quick Peek at the Data\n","We first visualize all the data we have to get a rough idea about how the data looks like.\n","\n","As we can see in the plot below, the hourly temperatures follow a clear seasonal trend. Wind speed, pressure, humidity and precipitation data seem to have much higher variance and randomness.\n","\n","It might be more meaningful to make a model to predict temperature, rather than some of the other more noisy data columns.\n"]},{"cell_type":"code","id":"891fa12a-11e5-453e-aaef-c1337775960a","metadata":{},"outputs":[],"source":["# Quick overview of columns\nplt.figure(figsize=(30, 12))\ni = 1\nfor col in plot_cols:\n plt.subplot(len(plot_cols), 1, i)\n plt.plot(data[col].values)\n plt.title(col)\n i += 1\nplt.subplots_adjust(hspace=0.5)\nplt.show()\n"]},{"cell_type":"markdown","id":"1c7db741-79d9-4f5b-b886-23f7c2bb11d7","metadata":{},"outputs":[],"source":["#### Feature Dependencies\n","Now we explore how the features (columns) of our data are related to each other. This helps in deciding which features to use when modelling a classifier or regresser. We ideally want independent features to be classified independently and likewise dependent features to be contributing to the same model.\n","\n","We can see from the correlation plots how some features are somewhat correlated and could be used as additional data (perhaps for augmenting) when training a classifier.\n"]},{"cell_type":"code","id":"5583332f-818a-4e1a-aad4-6d33db503892","metadata":{},"outputs":[],"source":["# Plot correlation matrix\nf, ax = plt.subplots(figsize=(7, 7))\ncorr = data[plot_cols].corr()\nsns.heatmap(corr, mask=np.zeros_like(corr, dtype=np.bool), cmap=sns.diverging_palette(220, 10, as_cmap=True), square=True, ax=ax);"]},{"cell_type":"markdown","id":"c1f92671-4bbe-42d5-ac1c-307a3b222c65","metadata":{},"outputs":[],"source":["Additionally we also visualize the joint distrubitions in the form of pairplots/scatter plots to see (qualitatively) the way in which these features are related in more detail over just the correlation. They are essentially 2D joint distributions in the case of off-diagonal subplots and the histogram (an approximation to the probability distribution) in case of the diagonal subplots.\n"]},{"cell_type":"code","id":"3ea317bf-69e6-4843-bba1-00bd8b490c21","metadata":{},"outputs":[],"source":["# Plot pairplots\nsns.pairplot(data[plot_cols])"]},{"cell_type":"markdown","id":"ad81d7bb-e509-4bba-97cd-0e2f26d6b269","metadata":{},"outputs":[],"source":["## Analyze Trends in the Data\n","Now that we have explored the whole dataset and the features on a high level, let us focus on one particular feature - dry_bulb_temp_f, the dry bulb temperature in degrees Fahrenheit. This is what we mean when we refer to \"air temperature\". This is the most common feature used in temperature prediction, and here we explore it in further detail.\n","\n","We first start with plotting the data for all 9 years in monthly buckets then drill down to a single year to notice (qualitatively) the overall trend in the data. We can see from the plots that every year has roughly a sinousoidal nature to the temperature with some anomalies around 2013-2014. Upon further drilling down we see that each year's data is not the smooth sinousoid but rather a jagged and noisy one. But the overall trend still is a sinousoid.\n"]},{"cell_type":"code","id":"f692697d-5526-4c36-9f0a-fa47640e6080","metadata":{},"outputs":[],"source":["plt.figure(figsize=(15,7))\n\nTEMP_COL = 'dry_bulb_temp_f'\n# Plot temperature data converted to a monthly frequency \nplt.subplot(1, 2, 1)\n#data = data[~data.index.duplicated(keep='first')]\ndata[TEMP_COL].asfreq('ME').plot() \nplt.title('Monthly Temperature')\nplt.ylabel('Temperature (F)')\n\n# Zoom in on a year and plot temperature data converted to a daily frequency \nplt.subplot(1, 2, 2)\ndata.loc['2017', TEMP_COL].asfreq('D').plot()\nplt.title('Daily Temperature in 2017')\nplt.ylabel('Temperature (F)')\n\nplt.tight_layout()\nplt.show()"]},{"cell_type":"markdown","id":"ae77d48c-4aa7-4e40-babd-34c0530d3a08","metadata":{},"outputs":[],"source":["Next, we plot the change (delta) in temperature and notice that it is lowest around the middle of the year. That is expected behaviour as the gradient of the sinousoid near it's peak is zero.\n"]},{"cell_type":"code","id":"cd5dcd71-fad5-4ad0-a18b-3e0d359bc8e3","metadata":{},"outputs":[],"source":["plt.figure(figsize=(15,7))\n\n# Plot percent change of daily temperature in 2017\nplt.subplot(1, 2, 1)\ndata.loc['2017', TEMP_COL].asfreq('D').div(\n data.loc['2017', TEMP_COL].asfreq('D').shift()\n).plot()\nplt.title('% Change in Daily Temperature in 2017')\nplt.ylabel('% Change')\n\n# Plot absolute change of temperature in 2017 with daily frequency\nplt.subplot(1, 2, 2)\n\ndata.loc['2017', TEMP_COL].asfreq('D').diff().plot()\n\nplt.title('Absolute Change in Daily Temperature in 2017')\nplt.ylabel('Temperature (F)')\n\nplt.tight_layout()\nplt.show()"]},{"cell_type":"markdown","id":"b4866b90-b0ff-4c89-b57b-d327550869c3","metadata":{},"outputs":[],"source":["Finally we apply some smoothing to the data in the form of a rolling/moving average. This is the simplest form of de-noising the data. As we can see from the plots, the average (plotted in blue) roughly traces the sinousoid and is now much smoother. This can improve the accuracy of a regression model trained to predict temperatures within a reasonable margin of error.\n"]},{"cell_type":"code","id":"66b2971d-c5da-4a06-a29d-7df65f9cea4a","metadata":{},"outputs":[],"source":["plt.figure(figsize=(15, 7))\n\n# Plot rolling mean of temperature \nplt.subplot(1, 2, 1)\ndata.loc['2017', TEMP_COL].rolling('5D').mean().plot(zorder=2) # 5-day rolling mean\ndata.loc['2017', TEMP_COL].plot(zorder=1)\nplt.legend(['Rolling', 'Temp'])\nplt.title('Rolling Avg in Hourly Temperature in 2017')\nplt.ylabel('Temperature (F)')\n\n# Plot rolling mean of temperature for Jan–Mar 2017\nplt.subplot(1, 2, 2)\ndata.loc['2017-01':'2017-03', TEMP_COL].rolling('2D').mean().plot(zorder=2) # 2-day rolling mean\ndata.loc['2017-01':'2017-03', TEMP_COL].plot(zorder=1)\nplt.legend(['Rolling', 'Temp'])\nplt.title('Rolling Avg in Hourly Temperature in Winter 2017')\nplt.ylabel('Temperature (F)')\n\nplt.tight_layout()\nplt.show()\n\n"]},{"cell_type":"markdown","id":"d4fb0f1b-0dc1-43ab-a9cd-b9f2762538a3","metadata":{},"outputs":[],"source":["## Part 3: Time Series Forecasting\n"]},{"cell_type":"code","id":"00899ddb-14cd-4b79-a72d-7394c3abf395","metadata":{},"outputs":[],"source":["%%capture\n# Install required libraries\n!pip install scikit-learn\n!!pip install statsmodels"]},{"cell_type":"code","id":"e84305b2-cbfb-40b8-a000-221620924f72","metadata":{},"outputs":[],"source":["#Import required modules\nfrom sklearn.metrics import mean_squared_error\nfrom statsmodels.tsa.statespace.sarimax import SARIMAX\nfrom matplotlib import pyplot as plt"]},{"cell_type":"markdown","id":"f659056e-c0dd-4065-aeed-40a56df4f0f6","metadata":{},"outputs":[],"source":["For purposes of time-series modeling, we will restrict our analysis to a 2-year sample of the dataset to avoid overly long model-training times.\n"]},{"cell_type":"code","id":"f7319fb2-ebb6-4fc8-9446-28d5627dad5b","metadata":{},"outputs":[],"source":["data.reset_index(inplace=True)\ndata['DATE'] = pd.to_datetime(data['DATE'])\n\n# Step 2: Set 'DATE' as the index\ndata.set_index('DATE', inplace=True)\n\n# Step 3: Now slice using date strings\nsample = data['2016-01-01':'2018-01-01']\nsample.info()"]},{"cell_type":"markdown","id":"d7558868-f04a-463c-bea0-08e1849b2d76","metadata":{},"outputs":[],"source":["### Create Training/Validation/Test Splits\n","Before we attempt any time-series analysis and prediction, we should split the dataset into training, validation and test sets. We use a portion of the data for training, and a portion of future data for our validation and test sets.\n","\n","If we instead trained a model on the full dataset, the model would learn to be very good at making predictions on that particular dataset, essentially just copying the answers it knows. However, when presented with data the model has not seen , it would perform poorly since it has not learned how to generalize its answers.\n","\n","By training on a portion of the dataset and testing the model's performance on another portion of the dataset (which data the model has not seen in training), we try to avoid our models \"over-fitting\" the dataset and make them better at predicting temperatures given unseen, future data. This process of splitting the dataset and evaluating a model's performance on the validation and test sets is commonly known as cross-validation).\n","\n","By default here we use 80% of the data for the training set and 10% each for validation and test sets.\n"]},{"cell_type":"code","id":"bf3f4a04-0970-47a5-a605-9cf27e9c9c70","metadata":{},"outputs":[],"source":["def split_data(data, val_size=0.1, test_size=0.1):\n \"\"\"\n Splits data to training, validation and testing parts\n \"\"\"\n ntest = int(round(len(data) * (1 - test_size)))\n nval = int(round(len(data) * (1 - test_size - val_size)))\n\n df_train, df_val, df_test = data.iloc[:nval], data.iloc[nval:ntest], data.iloc[ntest:]\n \n return df_train, df_val, df_test\n\n\n# Create data split\ndf_train, df_val, df_test = split_data(sample)\n\nprint('Total data size: {} rows'.format(len(sample)))\nprint('Training set size: {} rows'.format(len(df_train)))\nprint('Validation set size: {} rows'.format(len(df_val)))\nprint('Test set size: {} rows'.format(len(df_test)))"]},{"cell_type":"markdown","id":"6020874f-c208-4192-bf89-a1cd3526dbcc","metadata":{},"outputs":[],"source":["## Explore Baseline Models\n","In this section, we'll create a few simple predictive models of temperature, using shifting and rolling averages. These will serve as a baseline against which we can compare more sophisticated models.\n","\n","Using values at recent timesteps (such as the most recent timestep t-1 and second-most recent timestep t-2) to predict the current value at time t is what's known as persistence modeling, or using the last observed value to predict the next following value. These preceding timesteps are often referred to in time-series analysis as lags. So, the value at time t-1 is known as the 1st lag and the value at time t-2 is the 2nd lag.\n","\n","We can also create baselines based on rolling (or moving) averages. This is a time-series constructed by averaging each lagged value up to the selected lag. For example, a 6-period (or 6-lag) rolling avearge is the average of the previous 6 hourly lags t-6 to t-1.\n","\n","Our baseline models will be:\n","\n","- 1st lag - i.e. values at t-1\n","- 2nd lag - i.e. values at t-2\n","- 6-lag rolling average\n","- 12-lag rolling average\n"]},{"cell_type":"code","id":"a0eea160-f154-4734-bc71-a66f7ed11447","metadata":{},"outputs":[],"source":["%%capture\n# define the column containing the data we wish to model - in this case Dry Bulb Temperature (F)\nY_COL = 'dry_bulb_temp_f'\n\n# Use shifting and rolling averages to predict Y_COL (t)\nn_in = 2\nn_out = 1\nfeatures = [Y_COL]\nn_features = len(features)\n\n# create the baseline on the entire sample dataset.\n# we will evaluate the prediction error on the validation set\nbaseline = sample[[Y_COL]].loc[:]\nbaseline['{} (t-1)'.format(Y_COL)] = baseline[Y_COL].shift(1)\nbaseline['{} (t-2)'.format(Y_COL)] = baseline[Y_COL].shift(2)\nbaseline['{} (6hr rollavg)'.format(Y_COL)] = baseline[Y_COL].rolling('6H').mean()\nbaseline['{} (12hr rollavg)'.format(Y_COL)] = baseline[Y_COL].rolling('12H').mean()\nbaseline.dropna(inplace=True)\nbaseline.head(10)"]},{"cell_type":"markdown","id":"e2f5545d-ba36-4588-b67c-8b8d85bc6249","metadata":{},"outputs":[],"source":["Next, we will plot data from our validation dataset to get a sense for how well these baseline models predict the next hourly temperature. Note thatd we only use a few days of data in order to make the plot easier to view.\n"]},{"cell_type":"code","id":"09feead1-b1e9-48a4-931a-dfed8ca60638","metadata":{},"outputs":[],"source":["# plot first 7 days of the validation set, 168 hours \nstart = df_val.index[0]\nend = df_val.index[167]\nsliced = baseline[start:end]"]},{"cell_type":"code","id":"535e71c4-8473-4fbf-a91c-c0389d85a24f","metadata":{},"outputs":[],"source":["# Plot baseline predictions sample\ncols = ['dry_bulb_temp_f', 'dry_bulb_temp_f (t-1)', 'dry_bulb_temp_f (t-2)', 'dry_bulb_temp_f (6hr rollavg)', 'dry_bulb_temp_f (12hr rollavg)']\nsliced[cols].plot()\n\nplt.legend(['t', 't-1', 't-2', '6hr', '12hr'], loc=2, ncol=3)\nplt.title('Baselines for First 7 Days of Validation Set')\nplt.ylabel('Temperature (F)')\nplt.tight_layout()\nplt.rcParams['figure.dpi'] = 100\nplt.show()"]},{"cell_type":"markdown","id":"b0ccdd3f-fdd1-40c7-bf42-f96e1b4602cc","metadata":{},"outputs":[],"source":["Evaluate baseline models\n","As you can perhaps see from the graph above, the lagged baselines appear to do a better job of forecasting temperatures than the rolling average baselines.\n","\n","In order to evaluate our baseline models more precisely, we need to answer the question \"how well do our models predict future temperature?\". In regression problems involving prediction of a numerical value, we often use a measure of the difference between our predicted value and the actual value. This is referred to as an error measure or error metric. A common measure is the Mean Squared Error (MSE):\n","\n","$$\n","\\text{MSE} = \\frac{1}{n} \\sum_{i=1}^{n} (y_i - \\hat{y}_i)^2\n","$$\n","\n","This is the average of the squared differences between predicted values \n"," and actual values \n"," .\n","\n","Because the MSE is in \"units squared\" it can be difficult to interpet, hence the Root Mean Squared Error (RMSE) is often used:\n","\n","$$\n","\\text{RMSE} = \\sqrt{\\frac{1}{n} \\sum_{i=1}^{n} (y_i - \\hat{y}_i)^2}\n","$$\n"," \n","This is the square root of the MSE, and is in the same units as the values \n"," . We can compare the RMSE (and MSE) values for different models and say that the model that has the lower MSE is better at predicting temperatures, all things equal. Note that MSE and RMSE will grow large quickly if the differences between predicted and actual values are large. This may or may not be a desired quality of your error measure. In this case, it is probably a good thing, since a model that makes large mistakes in temperature prediction will be much less useful than one which makes small mistakes.\n","\n","Next, we calculate the RMSE measure for each of our baseline models, on the full validation set.\n"]},{"cell_type":"code","id":"17631269-1b76-4298-b2e3-208c21596958","metadata":{},"outputs":[],"source":["# Calculating baseline RMSE\nstart_val = df_val.index[0]\nend_val = df_val.index[-1]\nbaseline_val = baseline[start_val:end_val]\n\nbaseline_y = baseline_val[Y_COL]\nbaseline_t1 = baseline_val['dry_bulb_temp_f (t-1)']\nbaseline_t2 = baseline_val['dry_bulb_temp_f (t-2)']\nbaseline_avg6 = baseline_val['dry_bulb_temp_f (6hr rollavg)']\nbaseline_avg12 = baseline_val['dry_bulb_temp_f (12hr rollavg)']\n\nrmse_t1 = round(np.sqrt(mean_squared_error(baseline_y, baseline_t1)), 2)\nrmse_t2 = round(np.sqrt(mean_squared_error(baseline_y, baseline_t2)), 2)\nrmse_avg6 = round(np.sqrt(mean_squared_error(baseline_y, baseline_avg6)), 2)\nrmse_avg12 = round(np.sqrt(mean_squared_error(baseline_y, baseline_avg12)), 2)\n\nprint('Baseline t-1 RMSE: {0:.3f}'.format(rmse_t1))\nprint('Baseline t-2 RMSE: {0:.3f}'.format(rmse_t2))\nprint('Baseline 6hr rollavg RMSE: {0:.3f}'.format(rmse_avg6))\nprint('Baseline 12hr rollavg RMSE: {0:.3f}'.format(rmse_avg12))"]},{"cell_type":"markdown","id":"22fd5bc7-ecc2-48f9-8b31-634cde8a72a8","metadata":{},"outputs":[],"source":["The RMSE results confirm what we saw in the graph above. It is clear that the rolling average baselines perform poorly. In fact, the t-2 lagged baseline is also not very good. It appears that the best baseline model is to simply use the current hour's temperature to predict the next hour's temperature!\n","\n","Can we do better than this simple baseline using more sophisticated models?\n"]},{"cell_type":"markdown","id":"1c72d2cf-19fa-4c11-aa01-1d167a9798b1","metadata":{},"outputs":[],"source":["## Train Statistical Time-series Analysis Models\n","In the previous section, we saw that a simple lag-1 baseline model performed reasonably well at forecasting temperature for the next hourly time step. This is perhaps not too surprising, given what we know about hourly temperatures. Generally, the temperature in a given hour will be quite closely related to the temperature in the previous hour. This phenomenon is very common in time-series analysis and is known as autocorrelation - that is, the time series is correlated with previous values of itself. More precisely, the values at time t are correlated with lagged values (which could be t-1, t-2 and so on).\n","\n","Another thing we saw previously is the concept of moving averages. In this case the moving-average baseline was not that good at prediction. However it is common in many time-series for a moving average to capture some of the underlying structure and be useful for prediction.\n","\n","In order to make our model better at predicting temperature, ideally we would want to take these aspects into account. Fortunately, the statistical community has a long history of analyzing time series and has created many different forecasting models.\n","\n","Here, we will explore one called SARIMAX - the Seasonal AutoRegressive Integrated Moving Average with eXogenous regressors model.\n","\n","This sounds like a very complex name, but if we look at the components of the name, we see that it includes autocorrelation (this is what auto regressive means) and moving averages, which are the components mentioned above.\n","\n","The SARIMAX model also allows including a seasonal model component as well as handling exogenous variables, which are external to the time-series value itself. For example, for temperature prediction we may wish to take into account not just previous temperature values, but perhaps other weather features which may have an effect on temperature (such as humidity, rainfall, wind, and so on).\n","\n","For the purposes of this notebook, we will not explore modeling of seasonal components or exogenous variables.\n","\n","If we drop the \"S\" and \"X\" from the model, we are left with an ARIMA model (Auto-regressive Integrated Moving Average). This is a very commonly used model for time-series analysis and we will use it in this notebook by only specifying the relevant model components of the full SARIMAX model.\n","\n","### Replicating a baseline model\n","As a starting point, we will see how we can use SARIMAX to create a simple model that in fact replicates one of the baselines we created previously. Auto-regression, as we have seen, means using values from preceding time periods to predict the current value. Recall that one of our baseline models was the 1st lag or t-1 model. In time-series analysis this is referred to as an AR(1) model, meaning an Auto-Regressive model for lag 1.\n","\n","Technically, the AR(1) model is not exactly the same as our baseline model. A statistical time series model like SARIMAX learns a set of weights to apply to each component of the model. These weights are set so as to best fit the dataset. We can think of our baseline as setting the weight for the t-1 lag to be exactly 1. In practice, our time-series model will not have a weight of exactly 1 (though it will likely be very close to that), hence the predictions will be slightly different.\n","\n","Now, lets fit our model to the dataset. First, we will set up the model inputs by taking the temperature column of our dataframe. We do this for training and validation sets.\n"]},{"cell_type":"code","id":"8529cc86-a8ea-431b-8b8d-c006709d3c02","metadata":{},"outputs":[],"source":["X_train = df_train[Y_COL]\nX_val = df_val[Y_COL]\nX_both = np.hstack((X_train, X_val))"]},{"cell_type":"markdown","id":"8a29fc92-caab-41f4-84c3-310e4313a9cd","metadata":{},"outputs":[],"source":["Here we created a variable called X_both to cover both the training and validation data. This is required later when we forecast values for our SARIMAX model, in order to give the model access to all the datapoints for which it must create forecasts. Note that the forecasts themselves will only be based on the model weights learned from the training data (this is important for over-fitting as we have seen above)!\n","\n","The SARIMAX model takes an argument called order: this specifies the components of the model and itself has 3 parts: (p, d, q). p denotes the lags for the AR model and q denotes the lags for the MA model. We will not cover the d parameter here. Taken together this specifies the parameters of the ARIMA model portion of SARIMAX.\n","\n","To create an AR(1) model, we set the order to be (1, 0, 0). This sets up the AR model to be a lag 1 model. Then, we fit our model on the training data and inspect a summary of the trained model.\n"]},{"cell_type":"code","id":"8b6f3020-ff94-4755-bf80-93dd69479f2a","metadata":{},"outputs":[],"source":["%%capture\norder = (1, 0, 0)\nmodel_ar1 = SARIMAX(X_train, order=order)\nresults_ar1 = model_ar1.fit()\nresults_ar1.summary()"]},{"cell_type":"markdown","id":"d6bb4e70-11bb-456a-b844-9e52327a57e8","metadata":{},"outputs":[],"source":["### Warnings:\n","[1] Covariance matrix calculated using the outer product of gradients (complex-step).\n","There's quite a lot of information printed out in the model summary above. Much of it is related to the statistical properties of our model.\n","\n","The most important thing for now is to look at the second table, where we can see a coef value of 0.9996 for the weight ar.L1. This tells us the model has set a weight for the 1st lag component of the AR model to be 0.9996. This is almost 1 and hence we should expect the prediction results to indeed be close to our t-1 baseline.\n","\n","Let's create our model forecast on the validation dataset. We will then plot a few data points like we did with our baseline models (using 7 days of validation data) and compute the RMSE value based on the full validation set.\n"]},{"cell_type":"code","id":"add39958-f5bd-4458-9bc8-a9cbfe44e374","metadata":{},"outputs":[],"source":["full_data_ar1 = SARIMAX(X_both, order=order)\nmodel_forecast_ar1 = full_data_ar1.filter(results_ar1.params)"]},{"cell_type":"code","id":"d0b215a1-d0c1-4f97-9a0a-f6419c705afc","metadata":{},"outputs":[],"source":["start = len(X_train)\nend = len(X_both)\nforecast_ar1 = model_forecast_ar1.predict(start=start, end=end - 1, dynamic=False)\n\n# plot actual vs predicted values for the same 7-day window for easier viewing\nplt.plot(sliced[Y_COL].values)\nplt.plot(forecast_ar1[:168], color='r', linestyle='--')\nplt.legend(['t', 'AR(1)'], loc=2)\nplt.title('AR(1) Model Predictions for First 7 Days of Validation Set')\nplt.ylabel('Temperature (F)')\nplt.tight_layout()\nplt.show()"]},{"cell_type":"markdown","id":"991a99f3-5e29-4751-8590-81524754c93e","metadata":{},"outputs":[],"source":["We can see that the plot looks almost identical to the plot above, for the t and t-1 baseline values.\n","\n","Next, we compute the RMSE values.\n"]},{"cell_type":"code","id":"8301ab71-1927-40eb-9555-58e524057f2d","metadata":{},"outputs":[],"source":["# compute print RMSE values\nrmse_ar1 = np.sqrt(mean_squared_error(baseline_val[Y_COL], forecast_ar1))\nprint('AR(1) RMSE: {0:.3f}'.format(rmse_ar1))\nprint('Baseline t-1 RMSE: {0:.3f}'.format(rmse_t1))"]},{"cell_type":"markdown","id":"928e24b4-5155-4270-b101-ea02a9ca242a","metadata":{},"outputs":[],"source":["We can see that the RMSE values for the validation set also almost identical.\n","\n","## Create a more complex model\n","One of our baseline models was a lag 2 model, i.e. t-2. We saw that it performed a lot worse than the t-1 baseline. Intuitively, this makes sense, since we are throwing away a lot of information about the most recent lag t-1. However, the t-2 lag still provides some useful information. In fact, for temperature prediction it's likely that the last few hours can provide some value.\n","\n","Fortunately, our ARIMA model framework provides an easy way to incorporate further lag information. We can construct a model that includes both the t-1 and t-2 lags. This is an AR(2) model (meaning an auto-regressive model up to lag 2). We can specify this with the model order parameter p=2.\n"]},{"cell_type":"code","id":"33c1a49f-2c11-4b2b-82df-91de61f51a02","metadata":{},"outputs":[],"source":["%%capture\norder = (2, 0, 0)\nmodel_ar2 = SARIMAX(X_train, order=order)\nresults_ar2 = model_ar2.fit()\nresults_ar2.summary()"]},{"cell_type":"markdown","id":"2697f624-10fd-4a08-aaa5-8f4c43ca2368","metadata":{},"outputs":[],"source":["#### Warnings:\n","[1] Covariance matrix calculated using the outer product of gradients (complex-step).\n","\n","This time, the results table indicates a weight for variable ar.L1 and ar.L2. Note the values are now quite different from 1 (or 0.5 say, for a simple equally-weighted model). Next, we compute the RMSE on the validation set.\n"]},{"cell_type":"code","id":"ebc34b7d-26a9-4118-9cf8-dc7af516091c","metadata":{},"outputs":[],"source":["full_data_ar2 = SARIMAX(X_both, order=order)\nmodel_forecast_ar2 = full_data_ar2.filter(results_ar2.params)\n\nstart = len(X_train)\nend = len(X_both)\nforecast_ar2 = model_forecast_ar2.predict(start=start, end=end - 1, dynamic=False)\n\n# compute print RMSE values\nrmse_ar2 = np.sqrt(mean_squared_error(baseline_val[Y_COL], forecast_ar2))\nprint('AR(2) RMSE: {0:.3f}'.format(rmse_ar2))\nprint('AR(1) RMSE: {0:.3f}'.format(rmse_ar1))\nprint('Baseline t-1 RMSE: {0:.3f}'.format(rmse_t1))"]},{"cell_type":"markdown","id":"5db6903c-4434-4b35-8a3a-2802db2da2ba","metadata":{},"outputs":[],"source":["We've improved the RMSE value by including information from the first two lags.\n","\n","In fact, you will see that if you continue to increase the p parameter value, the RMSE will continue to decrease, indicating that a few recent lags provide useful information to our model.\n","\n","### Incorporate moving averages\n","Finally, what if we also include moving average information in our model? The ARIMA framework makes this easy to do, by setting the order parameter q. A value of q=1 specifies a MA(1) model (including the first lag t-1), while q=6 would include all the lags from t-1 to t-6.\n","\n","Note that the moving average model component is a little different from the simple moving or rolling averages computed in the baseline models. The definition of the MA model is rather technical, but conceptually you can think of it as using a form of weighted moving average (compared to our baseline which would be a simple, unweighted average).\n","\n","Let's add an MA(1) component to our AR(2) model.\n"]},{"cell_type":"code","id":"b438e5df-8c84-4f4d-9dfc-5af9afe66e8b","metadata":{},"outputs":[],"source":["%%capture\norder = (2, 0, 1)\nmodel_ar2ma1 = SARIMAX(X_train, order=order)\nresults_ar2ma1 = model_ar2ma1.fit()\nresults_ar2ma1.summary()"]},{"cell_type":"markdown","id":"ad45c0d8-ab35-44ca-8a63-3d61022893cc","metadata":{},"outputs":[],"source":["#### Warnings:\n","[1] Covariance matrix calculated using the outer product of gradients (complex-step).\n","\n","We see the results table shows an additional weight value for ma.L1, our MA(1) component. Next, we compare the RMSE to the other models and finally plot all the model forecasts together - note we use a much smaller 48-hour window to make the plot readable for illustrative purposes.\n"]},{"cell_type":"code","id":"bf6b7835-7115-4838-ade0-c64a58e5d95e","metadata":{},"outputs":[],"source":["full_data_ar2ma1 = SARIMAX(X_both, order=order)\nmodel_forecast_ar2ma1 = full_data_ar2ma1.filter(results_ar2ma1.params)\n\nstart = len(X_train)\nend = len(X_both)\nforecast_ar2ma1 = model_forecast_ar2ma1.predict(start=start, end=end - 1, dynamic=False)\n\n# compute print RMSE values\nrmse_ar2ma1 = np.sqrt(mean_squared_error(baseline_val[Y_COL], forecast_ar2ma1))\nprint('AR(2) MA(1) RMSE: {0:.3f}'.format(rmse_ar2ma1))\nprint('AR(2) RMSE: {0:.3f}'.format(rmse_ar2))\nprint('AR(1) RMSE: {0:.3f}'.format(rmse_ar1))\nprint('Baseline t-1 RMSE: {0:.3f}'.format(rmse_t1))"]},{"cell_type":"code","id":"194ef45c-c0cd-49d7-9470-6e0a919e34fd","metadata":{},"outputs":[],"source":["# plot actual vs predicted values for a smaller 2-day window for easier viewing\nhrs = 48\nplt.plot(sliced[Y_COL][:hrs].values)\nplt.plot(forecast_ar1[:hrs], color='r', linestyle='--')\nplt.plot(forecast_ar2[:hrs], color='g', linestyle='--')\nplt.plot(forecast_ar2ma1[:hrs], color='c', linestyle='--')\nplt.legend(['t', 'AR(1)', 'AR(2)', 'AR(2) MA(1)'], loc=2, ncol=1)\nplt.title('ARIMA Model Predictions for First 48 hours of Validation Set')\nplt.ylabel('Temperature (F)')\nplt.tight_layout()\nplt.show()"]},{"cell_type":"markdown","id":"f5407d4e-a806-4590-a151-5930f134779a","metadata":{},"outputs":[],"source":["We've again managed to reduce the RMSE value for our model, indicating that adding the MA(1) component has improved our forecast!\n","\n","*Congratulations! You've applied the basics of time-series analysis for forecasting hourly temperatures. See if you can further improve the RMSE values by exploring the different values for the model parameters p, q and even d!*\n"]},{"cell_type":"markdown","id":"7a937f3f-31b1-4723-ad5f-4e87c659a657","metadata":{},"outputs":[],"source":["Copyright © 2019 IBM. \n"]}],"metadata":{"kernelspec":{"display_name":"Python 3 (ipykernel)","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.12.8"},"prev_pub_hash":"006aebe6050b9a272be12fc15a7d6ec732f1b0278054a70a5d6a12048547f79f"},"nbformat":4,"nbformat_minor":4}