{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.12.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[],"dockerImageVersionId":28755,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session\n\n# Use the kagglehub client library to attach Kaggle resources like competitions, datasets, and models to your session\n# Learn more about kagglehub: https://github.com/Kaggle/kagglehub/blob/main/README.md\n\nimport kagglehub\n# kagglehub.dataset_download('<owner>/<dataset-slug>')","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2026-09-17T08:54:43.082838Z","iopub.execute_input":"2026-09-17T08:54:43.083821Z","iopub.status.idle":"2026-09-17T08:54:44.726776Z","shell.execute_reply.started":"2026-09-17T08:54:43.083783Z","shell.execute_reply":"2026-09-17T08:54:44.725189Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"#importing libraries\nimport pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport warnings\n\nwarnings.filterwarnings(\"ignore\")\n\nprint(\"Libraries imported successfully!\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T08:54:44.729756Z","iopub.execute_input":"2026-09-17T08:54:44.73034Z","iopub.status.idle":"2026-09-17T08:54:45.901189Z","shell.execute_reply.started":"2026-09-17T08:54:44.730307Z","shell.execute_reply":"2026-09-17T08:54:45.899848Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"#creating datapath\n\n\ndata_path = \"/kaggle/input/competitions/favorita-grocery-sales-forecasting\"\n\nprint(os.listdir(data_path))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T08:54:55.080662Z","iopub.execute_input":"2026-09-17T08:54:55.081179Z","iopub.status.idle":"2026-09-17T08:54:55.090432Z","shell.execute_reply.started":"2026-09-17T08:54:55.081148Z","shell.execute_reply":"2026-09-17T08:54:55.089309Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"#check file size\nimport os\n\nfor file in os.listdir(data_path):\n    file_path = os.path.join(data_path, file)\n    size_mb = os.path.getsize(file_path) / (1024 ** 2)\n    print(f\"{file:30s} {size_mb:10.2f} MB\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T08:54:55.930179Z","iopub.execute_input":"2026-09-17T08:54:55.930975Z","iopub.status.idle":"2026-09-17T08:54:55.95054Z","shell.execute_reply.started":"2026-09-17T08:54:55.930942Z","shell.execute_reply":"2026-09-17T08:54:55.949526Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"#extracting data from the zips\n!7z e \"/kaggle/input/competitions/favorita-grocery-sales-forecasting/items.csv.7z\" -o\"/kaggle/working/favorita_data\" -y","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T08:54:56.412369Z","iopub.execute_input":"2026-09-17T08:54:56.413322Z","iopub.status.idle":"2026-09-17T08:54:56.584803Z","shell.execute_reply.started":"2026-09-17T08:54:56.413272Z","shell.execute_reply":"2026-09-17T08:54:56.583578Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"working_path = \"/kaggle/working/favorita_data\"\n\nprint(os.listdir(working_path))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T08:54:56.870362Z","iopub.execute_input":"2026-09-17T08:54:56.870744Z","iopub.status.idle":"2026-09-17T08:54:56.877804Z","shell.execute_reply.started":"2026-09-17T08:54:56.870706Z","shell.execute_reply":"2026-09-17T08:54:56.876529Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import pandas as pd\n\nitems = pd.read_csv(\n    \"/kaggle/working/favorita_data/items.csv\"\n)\n\nprint(\"Shape:\", items.shape)\n\ndisplay(items.head())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T08:55:01.116115Z","iopub.execute_input":"2026-09-17T08:55:01.117111Z","iopub.status.idle":"2026-09-17T08:55:01.176566Z","shell.execute_reply.started":"2026-09-17T08:55:01.117071Z","shell.execute_reply":"2026-09-17T08:55:01.175181Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!7z","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T08:55:01.329759Z","iopub.execute_input":"2026-09-17T08:55:01.330326Z","iopub.status.idle":"2026-09-17T08:55:01.464069Z","shell.execute_reply.started":"2026-09-17T08:55:01.330241Z","shell.execute_reply":"2026-09-17T08:55:01.462225Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"items.dtypes","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T08:55:02.310633Z","iopub.execute_input":"2026-09-17T08:55:02.311547Z","iopub.status.idle":"2026-09-17T08:55:02.322467Z","shell.execute_reply.started":"2026-09-17T08:55:02.311497Z","shell.execute_reply":"2026-09-17T08:55:02.321304Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"items.shape","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T08:55:05.67147Z","iopub.execute_input":"2026-09-17T08:55:05.672142Z","iopub.status.idle":"2026-09-17T08:55:05.68165Z","shell.execute_reply.started":"2026-09-17T08:55:05.672107Z","shell.execute_reply":"2026-09-17T08:55:05.679915Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"items.nunique()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T08:55:05.98766Z","iopub.execute_input":"2026-09-17T08:55:05.988005Z","iopub.status.idle":"2026-09-17T08:55:06.007058Z","shell.execute_reply.started":"2026-09-17T08:55:05.987975Z","shell.execute_reply":"2026-09-17T08:55:06.005356Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"items.isnull().sum()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T08:55:06.999686Z","iopub.execute_input":"2026-09-17T08:55:07.000283Z","iopub.status.idle":"2026-09-17T08:55:07.016156Z","shell.execute_reply.started":"2026-09-17T08:55:07.000223Z","shell.execute_reply":"2026-09-17T08:55:07.013588Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"items[\"family\"].nunique()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T08:55:07.477333Z","iopub.execute_input":"2026-09-17T08:55:07.477774Z","iopub.status.idle":"2026-09-17T08:55:07.485569Z","shell.execute_reply.started":"2026-09-17T08:55:07.477746Z","shell.execute_reply":"2026-09-17T08:55:07.48441Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"items[\"family\"].value_counts()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T08:55:07.985385Z","iopub.execute_input":"2026-09-17T08:55:07.985724Z","iopub.status.idle":"2026-09-17T08:55:07.997905Z","shell.execute_reply.started":"2026-09-17T08:55:07.985697Z","shell.execute_reply":"2026-09-17T08:55:07.996825Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"items[\"perishable\"].value_counts()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T08:55:15.370163Z","iopub.execute_input":"2026-09-17T08:55:15.371784Z","iopub.status.idle":"2026-09-17T08:55:15.384102Z","shell.execute_reply.started":"2026-09-17T08:55:15.37174Z","shell.execute_reply":"2026-09-17T08:55:15.382909Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"items[\"perishable\"].value_counts(normalize=True) * 100","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T08:55:15.844868Z","iopub.execute_input":"2026-09-17T08:55:15.845489Z","iopub.status.idle":"2026-09-17T08:55:15.857779Z","shell.execute_reply.started":"2026-09-17T08:55:15.845452Z","shell.execute_reply":"2026-09-17T08:55:15.856373Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plt.figure(figsize=(12, 6))\n\nitems[\"family\"].value_counts().plot(kind=\"bar\")\n\nplt.title(\"Number of Products by Product Family\")\nplt.xlabel(\"Product Family\")\nplt.ylabel(\"Number of Products\")\nplt.xticks(rotation=85)\n\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T08:55:16.257437Z","iopub.execute_input":"2026-09-17T08:55:16.257826Z","iopub.status.idle":"2026-09-17T08:55:16.856388Z","shell.execute_reply.started":"2026-09-17T08:55:16.257798Z","shell.execute_reply":"2026-09-17T08:55:16.855445Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"#check stores\n!7z e \"/kaggle/input/competitions/favorita-grocery-sales-forecasting/stores.csv.7z\" \\\n    -o\"/kaggle/working/favorita_data\" -y","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T08:55:21.191025Z","iopub.execute_input":"2026-09-17T08:55:21.19141Z","iopub.status.idle":"2026-09-17T08:55:21.33449Z","shell.execute_reply.started":"2026-09-17T08:55:21.191375Z","shell.execute_reply":"2026-09-17T08:55:21.332731Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"stores = pd.read_csv(\n    \"/kaggle/working/favorita_data/stores.csv\"\n)\n\nprint(\"Shape:\", stores.shape)\ndisplay(stores.head())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T08:55:21.509552Z","iopub.execute_input":"2026-09-17T08:55:21.511048Z","iopub.status.idle":"2026-09-17T08:55:21.53239Z","shell.execute_reply.started":"2026-09-17T08:55:21.510984Z","shell.execute_reply":"2026-09-17T08:55:21.531192Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"stores.info()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T08:55:21.631001Z","iopub.execute_input":"2026-09-17T08:55:21.631357Z","iopub.status.idle":"2026-09-17T08:55:21.653822Z","shell.execute_reply.started":"2026-09-17T08:55:21.631329Z","shell.execute_reply":"2026-09-17T08:55:21.6528Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"stores.isnull().sum()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T08:55:21.786853Z","iopub.execute_input":"2026-09-17T08:55:21.78723Z","iopub.status.idle":"2026-09-17T08:55:21.799547Z","shell.execute_reply.started":"2026-09-17T08:55:21.787202Z","shell.execute_reply":"2026-09-17T08:55:21.798174Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"perishable_pct = (\n    items[\"perishable\"]\n    .value_counts(normalize=True)\n    .sort_index() * 100\n)\n\nprint(perishable_pct)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T08:55:21.976471Z","iopub.execute_input":"2026-09-17T08:55:21.976994Z","iopub.status.idle":"2026-09-17T08:55:21.993066Z","shell.execute_reply.started":"2026-09-17T08:55:21.976954Z","shell.execute_reply":"2026-09-17T08:55:21.991724Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!7z e \"/kaggle/input/competitions/favorita-grocery-sales-forecasting/stores.csv.7z\" \\\n    -o\"/kaggle/working/favorita_data\" -y","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T09:01:29.709031Z","iopub.execute_input":"2026-09-17T09:01:29.710269Z","iopub.status.idle":"2026-09-17T09:01:29.859789Z","shell.execute_reply.started":"2026-09-17T09:01:29.710196Z","shell.execute_reply":"2026-09-17T09:01:29.858289Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"stores = pd.read_csv(\n    \"/kaggle/working/favorita_data/stores.csv\"\n)\n\nprint(\"Shape:\", stores.shape)\n\ndisplay(stores.head())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T09:01:38.730053Z","iopub.execute_input":"2026-09-17T09:01:38.730496Z","iopub.status.idle":"2026-09-17T09:01:38.755605Z","shell.execute_reply.started":"2026-09-17T09:01:38.730454Z","shell.execute_reply":"2026-09-17T09:01:38.754792Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!7z e \"/kaggle/input/competitions/favorita-grocery-sales-forecasting/stores.csv.7z\" \\\n    -o\"/kaggle/working/favorita_data\" -y\n\nstores = pd.read_csv(\n    \"/kaggle/working/favorita_data/stores.csv\"\n)\n\nprint(\"Shape:\", stores.shape)\n\nprint(\"\\nFirst 5 rows:\")\ndisplay(stores.head())\n\nprint(\"\\nData types:\")\nprint(stores.dtypes)\n\nprint(\"\\nMissing values:\")\nprint(stores.isnull().sum())\n\nprint(\"\\nUnique values:\")\nprint(stores.nunique())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T09:01:50.740792Z","iopub.execute_input":"2026-09-17T09:01:50.741265Z","iopub.status.idle":"2026-09-17T09:01:50.901187Z","shell.execute_reply.started":"2026-09-17T09:01:50.741212Z","shell.execute_reply":"2026-09-17T09:01:50.899886Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plt.figure(figsize=(8, 5))\n\nstores[\"type\"].value_counts().plot(kind=\"bar\")\n\nplt.title(\"Number of Stores by Store Type\")\nplt.xlabel(\"Store Type\")\nplt.ylabel(\"Number of Stores\")\n\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T09:28:27.207292Z","iopub.execute_input":"2026-09-17T09:28:27.208266Z","iopub.status.idle":"2026-09-17T09:28:27.37867Z","shell.execute_reply.started":"2026-09-17T09:28:27.208201Z","shell.execute_reply":"2026-09-17T09:28:27.377419Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plt.figure(figsize=(10, 5))\n\nstores[\"city\"].value_counts().head(15).plot(kind=\"bar\")\n\nplt.title(\"Number of Stores by City\")\nplt.xlabel(\"City\")\nplt.ylabel(\"Number of Stores\")\n\nplt.xticks(rotation=45)\n\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T09:28:35.814094Z","iopub.execute_input":"2026-09-17T09:28:35.815037Z","iopub.status.idle":"2026-09-17T09:28:36.037584Z","shell.execute_reply.started":"2026-09-17T09:28:35.814994Z","shell.execute_reply":"2026-09-17T09:28:36.036361Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!7z l \"/kaggle/input/competitions/favorita-grocery-sales-forecasting/train.csv.7z\"","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T09:28:48.043928Z","iopub.execute_input":"2026-09-17T09:28:48.044351Z","iopub.status.idle":"2026-09-17T09:28:48.202465Z","shell.execute_reply.started":"2026-09-17T09:28:48.04432Z","shell.execute_reply":"2026-09-17T09:28:48.201102Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!df -h /kaggle/working","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T09:30:16.477458Z","iopub.execute_input":"2026-09-17T09:30:16.477883Z","iopub.status.idle":"2026-09-17T09:30:16.614583Z","shell.execute_reply.started":"2026-09-17T09:30:16.477843Z","shell.execute_reply":"2026-09-17T09:30:16.613025Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!7z e \"/kaggle/input/competitions/favorita-grocery-sales-forecasting/train.csv.7z\" \\\n    -o\"/kaggle/working/favorita_data\" -y","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T09:31:40.344721Z","iopub.execute_input":"2026-09-17T09:31:40.345394Z","iopub.status.idle":"2026-09-17T09:32:56.785226Z","shell.execute_reply.started":"2026-09-17T09:31:40.345347Z","shell.execute_reply":"2026-09-17T09:32:56.783692Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\n\ntrain_path = \"/kaggle/working/favorita_data/train.csv\"\n\nprint(os.path.exists(train_path))\n\nif os.path.exists(train_path):\n    size_gb = os.path.getsize(train_path) / (1024 ** 3)\n    print(f\"Train file size: {size_gb:.2f} GB\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T09:32:56.788216Z","iopub.execute_input":"2026-09-17T09:32:56.788653Z","iopub.status.idle":"2026-09-17T09:32:56.799142Z","shell.execute_reply.started":"2026-09-17T09:32:56.788613Z","shell.execute_reply":"2026-09-17T09:32:56.797352Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train_sample = pd.read_csv(\n    train_path,\n    nrows=10\n)\n\ndisplay(train_sample)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T09:33:25.347887Z","iopub.execute_input":"2026-09-17T09:33:25.348846Z","iopub.status.idle":"2026-09-17T09:33:25.377695Z","shell.execute_reply.started":"2026-09-17T09:33:25.348808Z","shell.execute_reply":"2026-09-17T09:33:25.376386Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"items[items[\"item_nbr\"] == 103665]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T09:34:49.652138Z","iopub.execute_input":"2026-09-17T09:34:49.658345Z","iopub.status.idle":"2026-09-17T09:34:49.689381Z","shell.execute_reply.started":"2026-09-17T09:34:49.658177Z","shell.execute_reply":"2026-09-17T09:34:49.687821Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train_sample.dtypes","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T09:35:03.529471Z","iopub.execute_input":"2026-09-17T09:35:03.530964Z","iopub.status.idle":"2026-09-17T09:35:03.54257Z","shell.execute_reply.started":"2026-09-17T09:35:03.530918Z","shell.execute_reply":"2026-09-17T09:35:03.540798Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train_sample['date'] = pd.to_datetime(train_sample['date'])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T09:39:11.671148Z","iopub.execute_input":"2026-09-17T09:39:11.673868Z","iopub.status.idle":"2026-09-17T09:39:11.734953Z","shell.execute_reply.started":"2026-09-17T09:39:11.673749Z","shell.execute_reply":"2026-09-17T09:39:11.733338Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train_sample.info()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T09:39:23.798552Z","iopub.execute_input":"2026-09-17T09:39:23.798987Z","iopub.status.idle":"2026-09-17T09:39:23.837144Z","shell.execute_reply.started":"2026-09-17T09:39:23.798956Z","shell.execute_reply":"2026-09-17T09:39:23.835574Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(train_sample[\"unit_sales\"].describe())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T09:39:52.108035Z","iopub.execute_input":"2026-09-17T09:39:52.109634Z","iopub.status.idle":"2026-09-17T09:39:52.134887Z","shell.execute_reply.started":"2026-09-17T09:39:52.109587Z","shell.execute_reply":"2026-09-17T09:39:52.13352Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train_path = \"/kaggle/working/favorita_data/train.csv\"\n\ntrain_sample = pd.read_csv(\n    train_path,\n    nrows=10\n)\n\ndisplay(train_sample)\n\nprint(\"\\nData types:\")\nprint(train_sample.dtypes)\n\nprint(\"\\nSummary:\")\nprint(train_sample[\"unit_sales\"].describe())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T09:40:13.326747Z","iopub.execute_input":"2026-09-17T09:40:13.327393Z","iopub.status.idle":"2026-09-17T09:40:13.383746Z","shell.execute_reply.started":"2026-09-17T09:40:13.327359Z","shell.execute_reply":"2026-09-17T09:40:13.382366Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train_sample[\"onpromotion\"].unique()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T09:42:12.496081Z","iopub.execute_input":"2026-09-17T09:42:12.496802Z","iopub.status.idle":"2026-09-17T09:42:12.510808Z","shell.execute_reply.started":"2026-09-17T09:42:12.49676Z","shell.execute_reply":"2026-09-17T09:42:12.50932Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train_sample.isnull().sum()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T09:42:38.746471Z","iopub.execute_input":"2026-09-17T09:42:38.746977Z","iopub.status.idle":"2026-09-17T09:42:38.756624Z","shell.execute_reply.started":"2026-09-17T09:42:38.746947Z","shell.execute_reply":"2026-09-17T09:42:38.755413Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train_sample = pd.read_csv(\n    train_path,\n    nrows=10000\n)\n\ntrain_sample[\"date\"] = pd.to_datetime(train_sample[\"date\"])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T09:43:17.278722Z","iopub.execute_input":"2026-09-17T09:43:17.279389Z","iopub.status.idle":"2026-09-17T09:43:17.313425Z","shell.execute_reply.started":"2026-09-17T09:43:17.279353Z","shell.execute_reply":"2026-09-17T09:43:17.312226Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(train_sample[\"onpromotion\"].unique())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T09:43:22.577183Z","iopub.execute_input":"2026-09-17T09:43:22.578761Z","iopub.status.idle":"2026-09-17T09:43:22.585771Z","shell.execute_reply.started":"2026-09-17T09:43:22.578703Z","shell.execute_reply":"2026-09-17T09:43:22.58431Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(train_sample[\"onpromotion\"].value_counts(dropna=False))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T09:43:30.455163Z","iopub.execute_input":"2026-09-17T09:43:30.455945Z","iopub.status.idle":"2026-09-17T09:43:30.465091Z","shell.execute_reply.started":"2026-09-17T09:43:30.45591Z","shell.execute_reply":"2026-09-17T09:43:30.464164Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(\"Minimum date:\", train_sample[\"date\"].min())\nprint(\"Maximum date:\", train_sample[\"date\"].max())\nprint(\"Number of unique dates:\", train_sample[\"date\"].nunique())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T09:43:55.791178Z","iopub.execute_input":"2026-09-17T09:43:55.791631Z","iopub.status.idle":"2026-09-17T09:43:55.810534Z","shell.execute_reply.started":"2026-09-17T09:43:55.791596Z","shell.execute_reply":"2026-09-17T09:43:55.809147Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(\"Number of stores:\", train_sample[\"store_nbr\"].nunique())\nprint(\"Number of products:\", train_sample[\"item_nbr\"].nunique())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T09:44:12.197955Z","iopub.execute_input":"2026-09-17T09:44:12.198282Z","iopub.status.idle":"2026-09-17T09:44:12.206189Z","shell.execute_reply.started":"2026-09-17T09:44:12.198239Z","shell.execute_reply":"2026-09-17T09:44:12.205084Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train_reader = pd.read_csv(\n    train_path,\n    usecols=[\"date\", \"onpromotion\"],\n    chunksize=100000\n)\n\nfor i, chunk in enumerate(train_reader):\n    \n    print(f\"\\nChunk {i}\")\n    print(\n        \"Date range:\",\n        chunk[\"date\"].min(),\n        \"to\",\n        chunk[\"date\"].max()\n    )\n    \n    print(\"Promotion values:\")\n    print(\n        chunk[\"onpromotion\"]\n        .value_counts(dropna=False)\n        .to_dict()\n    )\n    \n    if i == 4:\n        break","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T09:48:50.93844Z","iopub.execute_input":"2026-09-17T09:48:50.939679Z","iopub.status.idle":"2026-09-17T09:48:51.284675Z","shell.execute_reply.started":"2026-09-17T09:48:50.93964Z","shell.execute_reply":"2026-09-17T09:48:51.283326Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train_reader = pd.read_csv(\n    train_path,\n    usecols=[\"date\", \"onpromotion\"],\n    chunksize=100000\n)\n\nfor i, chunk in enumerate(train_reader):\n\n    if not chunk[\"onpromotion\"].isna().all():\n\n        print(\"First chunk containing promotion information:\")\n        print(\"Chunk:\", i)\n        print(\"Date range:\", chunk[\"date\"].min(), \"to\", chunk[\"date\"].max())\n\n        print(\"\\nPromotion values:\")\n        print(chunk[\"onpromotion\"].value_counts(dropna=False))\n\n        break","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T09:57:12.904468Z","iopub.execute_input":"2026-09-17T09:57:12.907594Z","iopub.status.idle":"2026-09-17T09:57:24.233191Z","shell.execute_reply.started":"2026-09-17T09:57:12.907544Z","shell.execute_reply":"2026-09-17T09:57:24.232063Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train_reader = pd.read_csv(\n    train_path,\n    usecols=[\"date\"],\n    chunksize=100000\n)\n\nmin_date = None\nmax_date = None\n\nfor chunk in train_reader:\n\n    chunk_min = chunk[\"date\"].min()\n    chunk_max = chunk[\"date\"].max()\n\n    if min_date is None or chunk_min < min_date:\n        min_date = chunk_min\n\n    if max_date is None or chunk_max > max_date:\n        max_date = chunk_max\n\nprint(\"Earliest date:\", min_date)\nprint(\"Latest date:\", max_date)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T09:58:30.62312Z","iopub.execute_input":"2026-09-17T09:58:30.623575Z","iopub.status.idle":"2026-09-17T09:59:53.970643Z","shell.execute_reply.started":"2026-09-17T09:58:30.623541Z","shell.execute_reply":"2026-09-17T09:59:53.969494Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"(pd.to_datetime(max_date) - pd.to_datetime(min_date)).days","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:02:07.954938Z","iopub.execute_input":"2026-09-17T10:02:07.955238Z","iopub.status.idle":"2026-09-17T10:02:07.966641Z","shell.execute_reply.started":"2026-09-17T10:02:07.955216Z","shell.execute_reply":"2026-09-17T10:02:07.965306Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train_reader = pd.read_csv(\n    train_path,\n    usecols=[\"date\"],\n    chunksize=100000\n)\n\ndates = set()\n\nfor chunk in train_reader:\n    dates.update(chunk[\"date\"].unique())\n\nprint(\"Number of unique dates:\", len(dates))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:02:41.664928Z","iopub.execute_input":"2026-09-17T10:02:41.665546Z","iopub.status.idle":"2026-09-17T10:03:52.729816Z","shell.execute_reply.started":"2026-09-17T10:02:41.66551Z","shell.execute_reply":"2026-09-17T10:03:52.728419Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"daily_sales_parts = []\n\ntrain_reader = pd.read_csv(\n    train_path,\n    usecols=[\"date\", \"unit_sales\"],\n    chunksize=100000\n)\n\nfor chunk in train_reader:\n    \n    daily_sales = (\n        chunk\n        .groupby(\"date\")[\"unit_sales\"]\n        .sum()\n    )\n    \n    daily_sales_parts.append(daily_sales)\n\ndaily_sales = (\n    pd.concat(daily_sales_parts)\n    .groupby(level=0)\n    .sum()\n    .reset_index()\n)\n\ndaily_sales.columns = [\"date\", \"total_sales\"]\n\ndaily_sales[\"date\"] = pd.to_datetime(daily_sales[\"date\"])\n\ndaily_sales = daily_sales.sort_values(\"date\")\n\ndisplay(daily_sales.head())\ndisplay(daily_sales.tail())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:04:15.658412Z","iopub.execute_input":"2026-09-17T10:04:15.658731Z","iopub.status.idle":"2026-09-17T10:05:30.523783Z","shell.execute_reply.started":"2026-09-17T10:04:15.658706Z","shell.execute_reply":"2026-09-17T10:05:30.522713Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(\"Shape:\", daily_sales.shape)\n\nprint(\"\\nDate range:\")\nprint(daily_sales[\"date\"].min(), \"to\", daily_sales[\"date\"].max())\n\nprint(\"\\nMissing values:\")\nprint(daily_sales.isnull().sum())\n\nprint(\"\\nSummary statistics:\")\nprint(daily_sales[\"total_sales\"].describe())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:05:30.525703Z","iopub.execute_input":"2026-09-17T10:05:30.526118Z","iopub.status.idle":"2026-09-17T10:05:30.54131Z","shell.execute_reply.started":"2026-09-17T10:05:30.526089Z","shell.execute_reply":"2026-09-17T10:05:30.540283Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plt.figure(figsize=(15, 6))\n\nplt.plot(\n    daily_sales[\"date\"],\n    daily_sales[\"total_sales\"]\n)\n\nplt.title(\"Daily Total Sales Over Time\")\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.xticks(rotation=45)\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:05:30.542634Z","iopub.execute_input":"2026-09-17T10:05:30.542998Z","iopub.status.idle":"2026-09-17T10:05:30.927199Z","shell.execute_reply.started":"2026-09-17T10:05:30.542971Z","shell.execute_reply":"2026-09-17T10:05:30.924167Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"full_dates = pd.date_range(\n    start=daily_sales[\"date\"].min(),\n    end=daily_sales[\"date\"].max(),\n    freq=\"D\"\n)\n\nmissing_dates = full_dates.difference(\n    daily_sales[\"date\"]\n)\n\nprint(\"Expected calendar days:\", len(full_dates))\nprint(\"Observed sales days:\", len(daily_sales))\nprint(\"Missing dates:\", len(missing_dates))\n\nprint(\"\\nMissing dates:\")\nprint(missing_dates)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:12:34.947189Z","iopub.execute_input":"2026-09-17T10:12:34.947865Z","iopub.status.idle":"2026-09-17T10:12:34.981954Z","shell.execute_reply.started":"2026-09-17T10:12:34.94783Z","shell.execute_reply":"2026-09-17T10:12:34.980671Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plt.figure(figsize=(15, 6))\n\nplt.plot(\n    daily_sales[\"date\"],\n    daily_sales[\"total_sales\"]\n)\n\nplt.title(\"Daily Total Sales Over Time\")\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.xticks(rotation=45)\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:13:12.069614Z","iopub.execute_input":"2026-09-17T10:13:12.070136Z","iopub.status.idle":"2026-09-17T10:13:12.511821Z","shell.execute_reply.started":"2026-09-17T10:13:12.070092Z","shell.execute_reply":"2026-09-17T10:13:12.510373Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\n# 1. Basic information\nprint(\"=\" * 60)\nprint(\"DAILY SALES DATASET\")\nprint(\"=\" * 60)\n\nprint(\"Shape:\", daily_sales.shape)\n\nprint(\"\\nDate range:\")\nprint(\n    daily_sales[\"date\"].min(),\n    \"to\",\n    daily_sales[\"date\"].max()\n)\n\nprint(\"\\nMissing values:\")\nprint(daily_sales.isnull().sum())\n\n\n# 2. Create the complete calendar range\nfull_dates = pd.date_range(\n    start=daily_sales[\"date\"].min(),\n    end=daily_sales[\"date\"].max(),\n    freq=\"D\"\n)\n\n# Find dates that are missing from our daily dataset\nmissing_dates = full_dates.difference(\n    daily_sales[\"date\"]\n)\n\nprint(\"\\n\" + \"=\" * 60)\nprint(\"MISSING CALENDAR DATES\")\nprint(\"=\" * 60)\n\nprint(\"Expected calendar days:\", len(full_dates))\nprint(\"Observed sales days:\", len(daily_sales))\nprint(\"Number of missing dates:\", len(missing_dates))\n\nprint(\"\\nMissing dates:\")\nprint(missing_dates)\n\n\n# 3. Summary statistics\nprint(\"\\n\" + \"=\" * 60)\nprint(\"SUMMARY STATISTICS\")\nprint(\"=\" * 60)\n\nprint(daily_sales[\"total_sales\"].describe())\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:14:38.394053Z","iopub.execute_input":"2026-09-17T10:14:38.395379Z","iopub.status.idle":"2026-09-17T10:14:38.417369Z","shell.execute_reply.started":"2026-09-17T10:14:38.395336Z","shell.execute_reply":"2026-09-17T10:14:38.416186Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\n# 4. Minimum-sales day\nmin_day = daily_sales.loc[\n    daily_sales[\"total_sales\"].idxmin()\n]\n\nprint(\"\\n\" + \"=\" * 60)\nprint(\"MINIMUM SALES DAY\")\nprint(\"=\" * 60)\n\nprint(min_day)\n\n\n# 5. Maximum-sales day\nmax_day = daily_sales.loc[\n    daily_sales[\"total_sales\"].idxmax()\n]\n\nprint(\"\\n\" + \"=\" * 60)\nprint(\"MAXIMUM SALES DAY\")\nprint(\"=\" * 60)\n\nprint(max_day)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:14:50.003666Z","iopub.execute_input":"2026-09-17T10:14:50.004441Z","iopub.status.idle":"2026-09-17T10:14:50.01858Z","shell.execute_reply.started":"2026-09-17T10:14:50.004402Z","shell.execute_reply":"2026-09-17T10:14:50.017204Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\n# 6. Plot daily total sales\nplt.figure(figsize=(15, 6))\n\nplt.plot(\n    daily_sales[\"date\"],\n    daily_sales[\"total_sales\"]\n)\n\nplt.title(\"Daily Total Sales Over Time\")\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.xticks(rotation=45)\nplt.tight_layout()\nplt.show()\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:14:57.817946Z","iopub.execute_input":"2026-09-17T10:14:57.818445Z","iopub.status.idle":"2026-09-17T10:14:58.135518Z","shell.execute_reply.started":"2026-09-17T10:14:57.818412Z","shell.execute_reply":"2026-09-17T10:14:58.134173Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\n# 7. Histogram of daily sales\nplt.figure(figsize=(10, 5))\n\nplt.hist(\n    daily_sales[\"total_sales\"],\n    bins=50\n)\n\nplt.title(\"Distribution of Daily Total Sales\")\nplt.xlabel(\"Total Unit Sales\")\nplt.ylabel(\"Frequency\")\n\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:15:08.457535Z","iopub.execute_input":"2026-09-17T10:15:08.457876Z","iopub.status.idle":"2026-09-17T10:15:08.69574Z","shell.execute_reply.started":"2026-09-17T10:15:08.45785Z","shell.execute_reply":"2026-09-17T10:15:08.694385Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# STEP 11 — WEEKLY AND MONTHLY SALES ANALYSIS\n# ============================================================\n\n# Make sure date is datetime\ndaily_sales[\"date\"] = pd.to_datetime(daily_sales[\"date\"])\n\n# Set date as index for resampling\ndaily_sales_indexed = daily_sales.set_index(\"date\")\n\n\n# ============================================================\n# 1. WEEKLY SALES\n# ============================================================\n\nweekly_sales = (\n    daily_sales_indexed[\"total_sales\"]\n    .resample(\"W\")\n    .sum()\n    .reset_index()\n)\n\nweekly_sales.columns = [\"date\", \"total_sales\"]\n\n\nprint(\"=\" * 60)\nprint(\"WEEKLY SALES\")\nprint(\"=\" * 60)\n\nprint(\"Shape:\", weekly_sales.shape)\n\nprint(\"\\nDate range:\")\nprint(\n    weekly_sales[\"date\"].min(),\n    \"to\",\n    weekly_sales[\"date\"].max()\n)\n\nprint(\"\\nSummary statistics:\")\nprint(weekly_sales[\"total_sales\"].describe())\n\n\n# ============================================================\n# 2. MONTHLY SALES\n# ============================================================\n\nmonthly_sales = (\n    daily_sales_indexed[\"total_sales\"]\n    .resample(\"ME\")\n    .sum()\n    .reset_index()\n)\n\nmonthly_sales.columns = [\"date\", \"total_sales\"]\n\n\nprint(\"\\n\" + \"=\" * 60)\nprint(\"MONTHLY SALES\")\nprint(\"=\" * 60)\n\nprint(\"Shape:\", monthly_sales.shape)\n\nprint(\"\\nDate range:\")\nprint(\n    monthly_sales[\"date\"].min(),\n    \"to\",\n    monthly_sales[\"date\"].max()\n)\n\nprint(\"\\nSummary statistics:\")\nprint(monthly_sales[\"total_sales\"].describe())\n\n\n# ============================================================\n# 3. WEEKLY SALES PLOT\n# ============================================================\n\nplt.figure(figsize=(15, 6))\n\nplt.plot(\n    weekly_sales[\"date\"],\n    weekly_sales[\"total_sales\"]\n)\n\nplt.title(\"Weekly Total Sales Over Time\")\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.xticks(rotation=45)\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 4. MONTHLY SALES PLOT\n# ============================================================\n\nplt.figure(figsize=(15, 6))\n\nplt.plot(\n    monthly_sales[\"date\"],\n    monthly_sales[\"total_sales\"]\n)\n\nplt.title(\"Monthly Total Sales Over Time\")\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.xticks(rotation=45)\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 5. AVERAGE SALES BY MONTH\n# ============================================================\n\nmonthly_sales[\"month\"] = monthly_sales[\"date\"].dt.month\n\navg_sales_by_month = (\n    monthly_sales\n    .groupby(\"month\")[\"total_sales\"]\n    .mean()\n)\n\nprint(\"\\n\" + \"=\" * 60)\nprint(\"AVERAGE SALES BY MONTH\")\nprint(\"=\" * 60)\n\nprint(avg_sales_by_month)\n\n\n# ============================================================\n# 6. AVERAGE SALES BY DAY OF WEEK\n# ============================================================\n\ndaily_sales[\"day_of_week\"] = daily_sales[\"date\"].dt.dayofweek\n\navg_sales_by_day = (\n    daily_sales\n    .groupby(\"day_of_week\")[\"total_sales\"]\n    .mean()\n)\n\nprint(\"\\n\" + \"=\" * 60)\nprint(\"AVERAGE SALES BY DAY OF WEEK\")\nprint(\"=\" * 60)\n\nprint(avg_sales_by_day)\n\n\n# ============================================================\n# 7. MONTHLY SEASONALITY BAR CHART\n# ============================================================\n\nplt.figure(figsize=(10, 5))\n\nplt.bar(\n    avg_sales_by_month.index,\n    avg_sales_by_month.values\n)\n\nplt.title(\"Average Sales by Month\")\nplt.xlabel(\"Month\")\nplt.ylabel(\"Average Total Sales\")\n\nplt.xticks(range(1, 13))\n\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 8. DAY-OF-WEEK PATTERN\n# ============================================================\n\nplt.figure(figsize=(10, 5))\n\nplt.bar(\n    avg_sales_by_day.index,\n    avg_sales_by_day.values\n)\n\nplt.title(\"Average Sales by Day of Week\")\nplt.xlabel(\"Day of Week\")\nplt.ylabel(\"Average Total Sales\")\n\nplt.xticks(\n    range(7),\n    [\"Monday\", \"Tuesday\", \"Wednesday\",\n     \"Thursday\", \"Friday\", \"Saturday\", \"Sunday\"],\n    rotation=45\n)\n\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:16:01.526191Z","iopub.execute_input":"2026-09-17T10:16:01.527493Z","iopub.status.idle":"2026-09-17T10:16:03.255653Z","shell.execute_reply.started":"2026-09-17T10:16:01.527444Z","shell.execute_reply":"2026-09-17T10:16:03.254355Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# STEP 12 — OUTLIERS AND DATA QUALITY CHECK\n# ============================================================\n\nprint(\"=\" * 70)\nprint(\"STEP 12: OUTLIERS AND DATA QUALITY\")\nprint(\"=\" * 70)\n\n\n# ============================================================\n# 1. IQR OUTLIER DETECTION ON DAILY SALES\n# ============================================================\n\nQ1 = daily_sales[\"total_sales\"].quantile(0.25)\nQ3 = daily_sales[\"total_sales\"].quantile(0.75)\n\nIQR = Q3 - Q1\n\nlower_bound = Q1 - 1.5 * IQR\nupper_bound = Q3 + 1.5 * IQR\n\ndaily_outliers = daily_sales[\n    (daily_sales[\"total_sales\"] < lower_bound) |\n    (daily_sales[\"total_sales\"] > upper_bound)\n].copy()\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"DAILY SALES IQR OUTLIERS\")\nprint(\"=\" * 70)\n\nprint(\"Q1:\", Q1)\nprint(\"Q3:\", Q3)\nprint(\"IQR:\", IQR)\nprint(\"Lower bound:\", lower_bound)\nprint(\"Upper bound:\", upper_bound)\n\nprint(\"\\nNumber of outlier days:\", len(daily_outliers))\n\nprint(\"\\nOutlier days:\")\ndisplay(\n    daily_outliers\n    .sort_values(\"total_sales\")\n    .head(20)\n)\n\n\n# ============================================================\n# 2. LOWEST AND HIGHEST SALES DAYS\n# ============================================================\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"LOWEST 10 SALES DAYS\")\nprint(\"=\" * 70)\n\ndisplay(\n    daily_sales\n    .nsmallest(10, \"total_sales\")\n)\n\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"HIGHEST 10 SALES DAYS\")\nprint(\"=\" * 70)\n\ndisplay(\n    daily_sales\n    .nlargest(10, \"total_sales\")\n)\n\n\n# ============================================================\n# 3. CHECK ZERO AND NEGATIVE DAILY SALES\n# ============================================================\n\nzero_sales_days = (\n    daily_sales[\"total_sales\"] == 0\n).sum()\n\nnegative_sales_days = (\n    daily_sales[\"total_sales\"] < 0\n).sum()\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"DAILY SALES QUALITY\")\nprint(\"=\" * 70)\n\nprint(\"Zero-sales days:\", zero_sales_days)\nprint(\"Negative-sales days:\", negative_sales_days)\n\n\n# ============================================================\n# 4. SCAN RAW TRAIN DATA FOR NEGATIVE UNIT SALES\n# ============================================================\n\nnegative_rows = 0\nzero_rows = 0\ntotal_rows = 0\n\nmax_unit_sales = None\nmin_unit_sales = None\n\ntrain_reader = pd.read_csv(\n    train_path,\n    usecols=[\"unit_sales\"],\n    chunksize=100000\n)\n\nfor chunk in train_reader:\n\n    total_rows += len(chunk)\n\n    negative_rows += (\n        chunk[\"unit_sales\"] < 0\n    ).sum()\n\n    zero_rows += (\n        chunk[\"unit_sales\"] == 0\n    ).sum()\n\n    chunk_max = chunk[\"unit_sales\"].max()\n    chunk_min = chunk[\"unit_sales\"].min()\n\n    if max_unit_sales is None or chunk_max > max_unit_sales:\n        max_unit_sales = chunk_max\n\n    if min_unit_sales is None or chunk_min < min_unit_sales:\n        min_unit_sales = chunk_min\n\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"RAW UNIT SALES QUALITY CHECK\")\nprint(\"=\" * 70)\n\nprint(\"Rows scanned:\", total_rows)\nprint(\"Negative unit_sales rows:\", negative_rows)\nprint(\"Zero unit_sales rows:\", zero_rows)\nprint(\"Minimum unit_sales:\", min_unit_sales)\nprint(\"Maximum unit_sales:\", max_unit_sales)\n\n\n# ============================================================\n# 5. OUTLIER VISUALIZATION\n# ============================================================\n\nplt.figure(figsize=(12, 5))\n\nplt.boxplot(\n    daily_sales[\"total_sales\"],\n    vert=False\n)\n\nplt.title(\"Boxplot of Daily Total Sales\")\nplt.xlabel(\"Total Unit Sales\")\n\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 6. DAILY SALES WITH OUTLIERS HIGHLIGHTED\n# ============================================================\n\nplt.figure(figsize=(15, 6))\n\nplt.plot(\n    daily_sales[\"date\"],\n    daily_sales[\"total_sales\"],\n    label=\"Daily Sales\"\n)\n\nplt.scatter(\n    daily_outliers[\"date\"],\n    daily_outliers[\"total_sales\"],\n    label=\"IQR Outliers\"\n)\n\nplt.title(\"Daily Sales with Potential Outliers\")\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\nplt.xticks(rotation=45)\n\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:16:43.883439Z","iopub.execute_input":"2026-09-17T10:16:43.884944Z","iopub.status.idle":"2026-09-17T10:17:46.27141Z","shell.execute_reply.started":"2026-09-17T10:16:43.8849Z","shell.execute_reply":"2026-09-17T10:17:46.270569Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# STEP 13 — CREATE TIME-BASED + HOLIDAY FEATURES\n# ============================================================\n\n# ------------------------------------------------------------\n# 1. Create a complete daily calendar\n# ------------------------------------------------------------\n\nfull_dates = pd.date_range(\n    start=daily_sales[\"date\"].min(),\n    end=daily_sales[\"date\"].max(),\n    freq=\"D\"\n)\n\ndaily_features = pd.DataFrame({\n    \"date\": full_dates\n})\n\n# Merge our observed daily sales onto the complete calendar\ndaily_features = daily_features.merge(\n    daily_sales[[\"date\", \"total_sales\"]],\n    on=\"date\",\n    how=\"left\"\n)\n\n# Flag dates where no sales record exists\ndaily_features[\"sales_missing\"] = (\n    daily_features[\"total_sales\"].isna().astype(int)\n)\n\n\n# ------------------------------------------------------------\n# 2. Create time-based features\n# ------------------------------------------------------------\n\ndaily_features[\"year\"] = daily_features[\"date\"].dt.year\ndaily_features[\"month\"] = daily_features[\"date\"].dt.month\ndaily_features[\"quarter\"] = daily_features[\"date\"].dt.quarter\ndaily_features[\"day\"] = daily_features[\"date\"].dt.day\ndaily_features[\"day_of_week\"] = daily_features[\"date\"].dt.dayofweek\n\n# Monday = 0 ... Sunday = 6\ndaily_features[\"day_name\"] = (\n    daily_features[\"date\"].dt.day_name()\n)\n\ndaily_features[\"is_weekend\"] = (\n    daily_features[\"day_of_week\"] >= 5\n).astype(int)\n\n\n# ------------------------------------------------------------\n# 3. Extract holidays_events.csv\n# ------------------------------------------------------------\n\nholidays_path = (\n    \"/kaggle/input/competitions/\"\n    \"favorita-grocery-sales-forecasting/holidays_events.csv.7z\"\n)\n\n!7z e \"{holidays_path}\" \\\n    -o\"/kaggle/working/favorita_data\" -y\n\n\n# ------------------------------------------------------------\n# 4. Load holiday data\n# ------------------------------------------------------------\n\nholidays = pd.read_csv(\n    \"/kaggle/working/favorita_data/holidays_events.csv\"\n)\n\nholidays[\"date\"] = pd.to_datetime(\n    holidays[\"date\"]\n)\n\nprint(\"=\" * 70)\nprint(\"HOLIDAY DATASET\")\nprint(\"=\" * 70)\n\nprint(\"Shape:\", holidays.shape)\n\nprint(\"\\nColumns:\")\nprint(holidays.columns.tolist())\n\nprint(\"\\nFirst 5 rows:\")\ndisplay(holidays.head())\n\nprint(\"\\nMissing values:\")\nprint(holidays.isnull().sum())\n\n\n# ------------------------------------------------------------\n# 5. Create one holiday record per date\n# ------------------------------------------------------------\n\nholiday_info = (\n    holidays\n    .groupby(\"date\")\n    .agg(\n        holiday_count=(\"description\", \"size\"),\n        holiday_types=(\n            \"type\",\n            lambda x: \", \".join(\n                sorted(set(x.astype(str)))\n            )\n        ),\n        holiday_names=(\n            \"description\",\n            lambda x: \", \".join(\n                sorted(set(x.astype(str)))\n            )\n        ),\n        holiday_locales=(\n            \"locale\",\n            lambda x: \", \".join(\n                sorted(set(x.astype(str)))\n            )\n        )\n    )\n    .reset_index()\n)\n\nholiday_info[\"is_holiday\"] = 1\n\n\n# ------------------------------------------------------------\n# 6. Merge holiday information with daily dataset\n# ------------------------------------------------------------\n\ndaily_features = daily_features.merge(\n    holiday_info,\n    on=\"date\",\n    how=\"left\"\n)\n\n\n# ------------------------------------------------------------\n# 7. Fill holiday indicators\n# ------------------------------------------------------------\n\ndaily_features[\"is_holiday\"] = (\n    daily_features[\"is_holiday\"]\n    .fillna(0)\n    .astype(int)\n)\n\ndaily_features[\"holiday_count\"] = (\n    daily_features[\"holiday_count\"]\n    .fillna(0)\n    .astype(int)\n)\n\ndaily_features[\"holiday_types\"] = (\n    daily_features[\"holiday_types\"]\n    .fillna(\"None\")\n)\n\ndaily_features[\"holiday_names\"] = (\n    daily_features[\"holiday_names\"]\n    .fillna(\"None\")\n)\n\ndaily_features[\"holiday_locales\"] = (\n    daily_features[\"holiday_locales\"]\n    .fillna(\"None\")\n)\n\n\n# ------------------------------------------------------------\n# 8. Display the final feature table\n# ------------------------------------------------------------\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"DAILY FEATURE DATASET\")\nprint(\"=\" * 70)\n\nprint(\"Shape:\", daily_features.shape)\n\nprint(\"\\nColumns:\")\nprint(daily_features.columns.tolist())\n\nprint(\"\\nFirst 10 rows:\")\ndisplay(daily_features.head(10))\n\nprint(\"\\nMissing sales dates:\")\ndisplay(\n    daily_features[\n        daily_features[\"sales_missing\"] == 1\n    ]\n)\n\nprint(\"\\nHoliday dates among missing-sales dates:\")\ndisplay(\n    daily_features[\n        (daily_features[\"sales_missing\"] == 1) &\n        (daily_features[\"is_holiday\"] == 1)\n    ]\n)\n\n\n# ------------------------------------------------------------\n# 9. Check feature distributions\n# ------------------------------------------------------------\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"TIME FEATURE SUMMARY\")\nprint(\"=\" * 70)\n\nprint(\"\\nYear counts:\")\nprint(daily_features[\"year\"].value_counts().sort_index())\n\nprint(\"\\nMonth counts:\")\nprint(daily_features[\"month\"].value_counts().sort_index())\n\nprint(\"\\nDay-of-week counts:\")\nprint(\n    daily_features[\"day_of_week\"]\n    .value_counts()\n    .sort_index()\n)\n\nprint(\"\\nWeekend counts:\")\nprint(\n    daily_features[\"is_weekend\"]\n    .value_counts()\n)\n\nprint(\"\\nHoliday counts:\")\nprint(\n    daily_features[\"is_holiday\"]\n    .value_counts()\n)\n\n\n# ------------------------------------------------------------\n# 10. Save the current feature dataset\n# ------------------------------------------------------------\n\ndaily_features.to_csv(\n    \"/kaggle/working/daily_sales_features.csv\",\n    index=False\n)\n\nprint(\"\\nSaved:\")\nprint(\"/kaggle/working/daily_sales_features.csv\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:18:25.256768Z","iopub.execute_input":"2026-09-17T10:18:25.257241Z","iopub.status.idle":"2026-09-17T10:18:25.673892Z","shell.execute_reply.started":"2026-09-17T10:18:25.257208Z","shell.execute_reply":"2026-09-17T10:18:25.672381Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# STEP 14 — TREND AND SEASONALITY ANALYSIS\n# ============================================================\n\n# Make sure dates are datetime\ndaily_features[\"date\"] = pd.to_datetime(daily_features[\"date\"])\n\n\n# ============================================================\n# 1. YEARLY TOTAL SALES\n# ============================================================\n\nyearly_sales = (\n    daily_features\n    .dropna(subset=[\"total_sales\"])\n    .groupby(\"year\")[\"total_sales\"]\n    .sum()\n    .reset_index()\n)\n\nprint(\"=\" * 70)\nprint(\"YEARLY TOTAL SALES\")\nprint(\"=\" * 70)\n\ndisplay(yearly_sales)\n\n\n# ============================================================\n# 2. YEARLY AVERAGE DAILY SALES\n# ============================================================\n\nyearly_avg_sales = (\n    daily_features\n    .dropna(subset=[\"total_sales\"])\n    .groupby(\"year\")[\"total_sales\"]\n    .mean()\n    .reset_index()\n)\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"YEARLY AVERAGE DAILY SALES\")\nprint(\"=\" * 70)\n\ndisplay(yearly_avg_sales)\n\n\n# ============================================================\n# 3. AVERAGE SALES BY MONTH\n# ============================================================\n\nmonthly_pattern = (\n    daily_features\n    .dropna(subset=[\"total_sales\"])\n    .groupby(\"month\")[\"total_sales\"]\n    .mean()\n)\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"AVERAGE DAILY SALES BY MONTH\")\nprint(\"=\" * 70)\n\ndisplay(monthly_pattern)\n\n\n# ============================================================\n# 4. AVERAGE SALES BY DAY OF WEEK\n# ============================================================\n\nweekday_pattern = (\n    daily_features\n    .dropna(subset=[\"total_sales\"])\n    .groupby(\n        [\"day_of_week\", \"day_name\"]\n    )[\"total_sales\"]\n    .mean()\n    .reset_index()\n)\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"AVERAGE DAILY SALES BY DAY OF WEEK\")\nprint(\"=\" * 70)\n\ndisplay(weekday_pattern)\n\n\n# ============================================================\n# 5. 30-DAY ROLLING AVERAGE\n# ============================================================\n\ndaily_features[\"rolling_30d\"] = (\n    daily_features[\"total_sales\"]\n    .rolling(window=30, min_periods=15)\n    .mean()\n)\n\n\n# ============================================================\n# 6. DAILY SALES + 30-DAY TREND\n# ============================================================\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    daily_features[\"date\"],\n    daily_features[\"total_sales\"],\n    alpha=0.35,\n    label=\"Daily Sales\"\n)\n\nplt.plot(\n    daily_features[\"date\"],\n    daily_features[\"rolling_30d\"],\n    linewidth=2,\n    label=\"30-Day Rolling Average\"\n)\n\nplt.title(\"Daily Sales with 30-Day Rolling Trend\")\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\nplt.xticks(rotation=45)\n\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 7. MONTHLY SEASONAL PATTERN\n# ============================================================\n\nplt.figure(figsize=(10, 5))\n\nplt.bar(\n    monthly_pattern.index,\n    monthly_pattern.values\n)\n\nplt.title(\"Average Daily Sales by Month\")\nplt.xlabel(\"Month\")\nplt.ylabel(\"Average Daily Sales\")\n\nplt.xticks(range(1, 13))\n\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 8. WEEKLY SEASONAL PATTERN\n# ============================================================\n\nplt.figure(figsize=(10, 5))\n\nplt.bar(\n    weekday_pattern[\"day_name\"],\n    weekday_pattern[\"total_sales\"]\n)\n\nplt.title(\"Average Daily Sales by Day of Week\")\nplt.xlabel(\"Day of Week\")\nplt.ylabel(\"Average Daily Sales\")\n\nplt.xticks(rotation=45)\n\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 9. YEAR × MONTH TABLE\n# ============================================================\n\nyear_month_sales = (\n    daily_features\n    .dropna(subset=[\"total_sales\"])\n    .groupby([\"year\", \"month\"])[\"total_sales\"]\n    .mean()\n    .unstack()\n)\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"YEAR × MONTH — AVERAGE DAILY SALES\")\nprint(\"=\" * 70)\n\ndisplay(year_month_sales)\n\n\n# ============================================================\n# 10. YEAR × MONTH HEATMAP\n# ============================================================\n\nplt.figure(figsize=(14, 6))\n\nsns.heatmap(\n    year_month_sales,\n    annot=True,\n    fmt=\".0f\"\n)\n\nplt.title(\"Average Daily Sales by Year and Month\")\nplt.xlabel(\"Month\")\nplt.ylabel(\"Year\")\n\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 11. YEARLY TREND PLOT\n# ============================================================\n\nplt.figure(figsize=(10, 5))\n\nplt.plot(\n    yearly_avg_sales[\"year\"],\n    yearly_avg_sales[\"total_sales\"],\n    marker=\"o\"\n)\n\nplt.title(\"Average Daily Sales by Year\")\nplt.xlabel(\"Year\")\nplt.ylabel(\"Average Daily Sales\")\n\nplt.xticks(\n    yearly_avg_sales[\"year\"]\n)\n\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 12. HOLIDAY VS NON-HOLIDAY SALES\n# ============================================================\n\nholiday_comparison = (\n    daily_features\n    .dropna(subset=[\"total_sales\"])\n    .groupby(\"is_holiday\")[\"total_sales\"]\n    .agg([\"count\", \"mean\", \"median\"])\n    .reset_index()\n)\n\nholiday_comparison[\"category\"] = (\n    holiday_comparison[\"is_holiday\"]\n    .map({\n        0: \"Non-Holiday\",\n        1: \"Holiday\"\n    })\n)\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"HOLIDAY VS NON-HOLIDAY SALES\")\nprint(\"=\" * 70)\n\ndisplay(holiday_comparison)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:19:34.430334Z","iopub.execute_input":"2026-09-17T10:19:34.430904Z","iopub.status.idle":"2026-09-17T10:19:35.858534Z","shell.execute_reply.started":"2026-09-17T10:19:34.430863Z","shell.execute_reply":"2026-09-17T10:19:35.857128Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# STEP 15 — FINALIZE MISSING VALUES + CLEAN FORECAST SERIES\n# ============================================================\n\n# Make a copy so we keep the original daily_features unchanged\nforecast_df = daily_features.copy()\n\n# ------------------------------------------------------------\n# 1. Preserve the original sales values\n# ------------------------------------------------------------\n\nforecast_df[\"original_total_sales\"] = forecast_df[\"total_sales\"]\n\n\n# ------------------------------------------------------------\n# 2. Identify structural missing sales days\n# ------------------------------------------------------------\n\nforecast_df[\"structural_missing\"] = (\n    forecast_df[\"total_sales\"].isna() &\n    (forecast_df[\"is_holiday\"] == 1)\n).astype(int)\n\n\n# ------------------------------------------------------------\n# 3. Fill structural missing sales with 0\n# ------------------------------------------------------------\n\nforecast_df[\"total_sales\"] = forecast_df[\"total_sales\"].fillna(0)\n\n\n# ------------------------------------------------------------\n# 4. Check whether any missing values remain\n# ------------------------------------------------------------\n\nprint(\"=\" * 70)\nprint(\"MISSING VALUE CHECK\")\nprint(\"=\" * 70)\n\nprint(\"Remaining missing total_sales:\",\n      forecast_df[\"total_sales\"].isna().sum())\n\nprint(\"Structural missing days:\",\n      forecast_df[\"structural_missing\"].sum())\n\nprint(\"\\nStructural missing dates:\")\ndisplay(\n    forecast_df[\n        forecast_df[\"structural_missing\"] == 1\n    ][[\n        \"date\",\n        \"original_total_sales\",\n        \"total_sales\",\n        \"holiday_names\",\n        \"is_holiday\",\n        \"structural_missing\"\n    ]]\n)\n\n\n# ------------------------------------------------------------\n# 5. Confirm the four Christmas dates\n# ------------------------------------------------------------\n\nchristmas_check = forecast_df[\n    (forecast_df[\"month\"] == 12) &\n    (forecast_df[\"day\"] == 25)\n][[\n    \"date\",\n    \"original_total_sales\",\n    \"total_sales\",\n    \"holiday_names\",\n    \"structural_missing\"\n]]\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"CHRISTMAS CHECK\")\nprint(\"=\" * 70)\n\ndisplay(christmas_check)\n\n\n# ------------------------------------------------------------\n# 6. IQR outlier flag\n# ------------------------------------------------------------\n\nQ1 = forecast_df[\"total_sales\"].quantile(0.25)\nQ3 = forecast_df[\"total_sales\"].quantile(0.75)\nIQR = Q3 - Q1\n\nlower_bound = Q1 - 1.5 * IQR\nupper_bound = Q3 + 1.5 * IQR\n\nforecast_df[\"outlier_flag\"] = (\n    (forecast_df[\"total_sales\"] < lower_bound) |\n    (forecast_df[\"total_sales\"] > upper_bound)\n).astype(int)\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"OUTLIER CHECK\")\nprint(\"=\" * 70)\n\nprint(\"Number of flagged outliers:\",\n      forecast_df[\"outlier_flag\"].sum())\n\nprint(\"\\nOutlier dates:\")\ndisplay(\n    forecast_df[\n        forecast_df[\"outlier_flag\"] == 1\n    ][[\n        \"date\",\n        \"total_sales\",\n        \"outlier_flag\"\n    ]]\n)\n\n\n# ------------------------------------------------------------\n# 7. Check final target series\n# ------------------------------------------------------------\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"FINAL FORECASTING SERIES\")\nprint(\"=\" * 70)\n\nprint(\"Shape:\", forecast_df.shape)\n\nprint(\"Date range:\",\n      forecast_df[\"date\"].min(),\n      \"to\",\n      forecast_df[\"date\"].max())\n\nprint(\"\\nMissing values:\")\nprint(\n    forecast_df[\n        [\"date\", \"total_sales\"]\n    ].isna().sum()\n)\n\nprint(\"\\nTarget summary:\")\nprint(\n    forecast_df[\"total_sales\"].describe()\n)\n\n\n# ------------------------------------------------------------\n# 8. Plot cleaned series\n# ------------------------------------------------------------\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    forecast_df[\"date\"],\n    forecast_df[\"total_sales\"]\n)\n\nplt.title(\"Cleaned Daily Total Sales\")\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.xticks(rotation=45)\nplt.tight_layout()\nplt.show()\n\n\n# ------------------------------------------------------------\n# 9. Save cleaned forecasting dataset\n# ------------------------------------------------------------\n\nforecast_df.to_csv(\n    \"/kaggle/working/final_forecasting_dataset.csv\",\n    index=False\n)\n\nprint(\"\\nSaved:\")\nprint(\"/kaggle/working/final_forecasting_dataset.csv\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:20:12.72121Z","iopub.execute_input":"2026-09-17T10:20:12.721641Z","iopub.status.idle":"2026-09-17T10:20:13.09674Z","shell.execute_reply.started":"2026-09-17T10:20:12.721606Z","shell.execute_reply":"2026-09-17T10:20:13.095604Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# STEP 16 — STATIONARITY TEST (ADF)\n# ============================================================\n\nfrom statsmodels.tsa.stattools import adfuller\n\n# ------------------------------------------------------------\n# 1. Prepare the time series\n# ------------------------------------------------------------\n\nts = (\n    forecast_df\n    .set_index(\"date\")[\"total_sales\"]\n    .asfreq(\"D\")\n)\n\nprint(\"=\" * 70)\nprint(\"TIME SERIES\")\nprint(\"=\" * 70)\n\nprint(\"Observations:\", len(ts))\nprint(\"Start:\", ts.index.min())\nprint(\"End:\", ts.index.max())\n\n\n# ------------------------------------------------------------\n# 2. Function for ADF test\n# ------------------------------------------------------------\n\ndef adf_test(series, name):\n    \n    result = adfuller(\n        series.dropna(),\n        autolag=\"AIC\"\n    )\n    \n    print(\"\\n\" + \"=\" * 70)\n    print(f\"ADF TEST — {name}\")\n    print(\"=\" * 70)\n    \n    print(\"ADF Statistic:\", result[0])\n    print(\"p-value:\", result[1])\n    print(\"Used Lags:\", result[2])\n    print(\"Number of Observations:\", result[3])\n    \n    print(\"\\nCritical Values:\")\n    \n    for key, value in result[4].items():\n        print(f\"   {key}: {value}\")\n    \n    if result[1] < 0.05:\n        print(\"\\nConclusion: Reject H0 → Series is stationary.\")\n    else:\n        print(\"\\nConclusion: Fail to reject H0 → Series is non-stationary.\")\n    \n    return result\n\n\n# ------------------------------------------------------------\n# 3. ADF on original series\n# ------------------------------------------------------------\n\nadf_original = adf_test(\n    ts,\n    \"Original Daily Sales\"\n)\n\n\n# ------------------------------------------------------------\n# 4. First-order differencing\n# ------------------------------------------------------------\n\nts_diff1 = ts.diff(1)\n\nadf_diff1 = adf_test(\n    ts_diff1,\n    \"First Difference\"\n)\n\n\n# ------------------------------------------------------------\n# 5. Seven-day seasonal difference\n# ------------------------------------------------------------\n\nts_diff7 = ts.diff(7)\n\nadf_diff7 = adf_test(\n    ts_diff7,\n    \"7-Day Seasonal Difference\"\n)\n\n\n# ------------------------------------------------------------\n# 6. Plot original vs first difference\n# ------------------------------------------------------------\n\nplt.figure(figsize=(15, 6))\n\nplt.plot(\n    ts.index,\n    ts.values\n)\n\nplt.title(\"Original Daily Sales\")\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Sales\")\n\nplt.tight_layout()\nplt.show()\n\n\nplt.figure(figsize=(15, 6))\n\nplt.plot(\n    ts_diff1.index,\n    ts_diff1.values\n)\n\nplt.axhline(\n    y=0,\n    linestyle=\"--\"\n)\n\nplt.title(\"First-Differenced Daily Sales\")\nplt.xlabel(\"Date\")\nplt.ylabel(\"Change in Sales\")\n\nplt.tight_layout()\nplt.show()\n\n\n# ------------------------------------------------------------\n# 7. Plot 7-day seasonal difference\n# ------------------------------------------------------------\n\nplt.figure(figsize=(15, 6))\n\nplt.plot(\n    ts_diff7.index,\n    ts_diff7.values\n)\n\nplt.axhline(\n    y=0,\n    linestyle=\"--\"\n)\n\nplt.title(\"7-Day Seasonal Difference\")\nplt.xlabel(\"Date\")\nplt.ylabel(\"7-Day Change in Sales\")\n\nplt.tight_layout()\nplt.show()\n\n\n# ------------------------------------------------------------\n# 8. Save the stationary/differenced series\n# ------------------------------------------------------------\n\nstationarity_df = pd.DataFrame({\n    \"date\": ts.index,\n    \"sales\": ts.values,\n    \"diff_1\": ts_diff1.values,\n    \"diff_7\": ts_diff7.values\n})\n\nstationarity_df.to_csv(\n    \"/kaggle/working/stationarity_analysis.csv\",\n    index=False\n)\n\nprint(\"\\nSaved:\")\nprint(\"/kaggle/working/stationarity_analysis.csv\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:20:46.161092Z","iopub.execute_input":"2026-09-17T10:20:46.162409Z","iopub.status.idle":"2026-09-17T10:20:47.72169Z","shell.execute_reply.started":"2026-09-17T10:20:46.162355Z","shell.execute_reply":"2026-09-17T10:20:47.720329Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# STEP 17 — ACF AND PACF ANALYSIS\n# ============================================================\n\nfrom statsmodels.graphics.tsaplots import plot_acf, plot_pacf\nfrom statsmodels.tsa.stattools import acf, pacf\n\n# ------------------------------------------------------------\n# 1. Prepare stationary series\n# ------------------------------------------------------------\n\n# First difference\ndiff1 = ts.diff(1).dropna()\n\n# 7-day seasonal difference\ndiff7 = ts.diff(7).dropna()\n\n\n# ============================================================\n# 2. ACF — FIRST DIFFERENCE\n# ============================================================\n\nplt.figure(figsize=(14, 5))\n\nplot_acf(\n    diff1,\n    lags=40,\n    ax=plt.gca()\n)\n\nplt.title(\"ACF — First-Differenced Daily Sales\")\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 3. PACF — FIRST DIFFERENCE\n# ============================================================\n\nplt.figure(figsize=(14, 5))\n\nplot_pacf(\n    diff1,\n    lags=40,\n    method=\"ywm\",\n    ax=plt.gca()\n)\n\nplt.title(\"PACF — First-Differenced Daily Sales\")\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 4. ACF — 7-DAY SEASONAL DIFFERENCE\n# ============================================================\n\nplt.figure(figsize=(14, 5))\n\nplot_acf(\n    diff7,\n    lags=60,\n    ax=plt.gca()\n)\n\nplt.title(\"ACF — 7-Day Seasonal Difference\")\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 5. PACF — 7-DAY SEASONAL DIFFERENCE\n# ============================================================\n\nplt.figure(figsize=(14, 5))\n\nplot_pacf(\n    diff7,\n    lags=60,\n    method=\"ywm\",\n    ax=plt.gca()\n)\n\nplt.title(\"PACF — 7-Day Seasonal Difference\")\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 6. NUMERICAL ACF VALUES\n# ============================================================\n\nacf_values = acf(\n    diff1,\n    nlags=14,\n    fft=True\n)\n\nprint(\"=\" * 70)\nprint(\"ACF VALUES — FIRST DIFFERENCE\")\nprint(\"=\" * 70)\n\nfor lag, value in enumerate(acf_values):\n    print(f\"Lag {lag:2d}: {value:.4f}\")\n\n\n# ============================================================\n# 7. NUMERICAL PACF VALUES\n# ============================================================\n\npacf_values = pacf(\n    diff1,\n    nlags=14,\n    method=\"ywm\"\n)\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"PACF VALUES — FIRST DIFFERENCE\")\nprint(\"=\" * 70)\n\nfor lag, value in enumerate(pacf_values):\n    print(f\"Lag {lag:2d}: {value:.4f}\")\n\n\n# ============================================================\n# 8. IMPORTANT SEASONAL LAGS\n# ============================================================\n\nacf_seasonal = acf(\n    diff7,\n    nlags=28,\n    fft=True\n)\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"IMPORTANT ACF LAGS — 7-DAY DIFFERENCE\")\nprint(\"=\" * 70)\n\nfor lag in [1, 2, 3, 6, 7, 8, 14, 21, 28]:\n    print(\n        f\"Lag {lag:2d}: \"\n        f\"{acf_seasonal[lag]:.4f}\"\n    )","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:21:13.955436Z","iopub.execute_input":"2026-09-17T10:21:13.956776Z","iopub.status.idle":"2026-09-17T10:21:14.800021Z","shell.execute_reply.started":"2026-09-17T10:21:13.956738Z","shell.execute_reply":"2026-09-17T10:21:14.799192Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# STEP 18 — CHRONOLOGICAL TRAIN / VALIDATION / TEST SPLIT\n# ============================================================\n\n# ------------------------------------------------------------\n# 1. Prepare the time series\n# ------------------------------------------------------------\n\nmodel_df = (\n    forecast_df[\n        [\"date\", \"total_sales\"]\n    ]\n    .copy()\n    .sort_values(\"date\")\n)\n\nmodel_df[\"date\"] = pd.to_datetime(\n    model_df[\"date\"]\n)\n\nmodel_df = model_df.set_index(\"date\")\n\n# Make sure we have a complete daily frequency\nmodel_df = model_df.asfreq(\"D\")\n\n\n# ------------------------------------------------------------\n# 2. Define split dates\n# ------------------------------------------------------------\n\ntrain_end = pd.Timestamp(\"2016-02-15\")\n\nvalidation_end = pd.Timestamp(\"2016-08-15\")\n\ntest_start = validation_end + pd.Timedelta(days=1)\n\n\n# ------------------------------------------------------------\n# 3. Create the three datasets\n# ------------------------------------------------------------\n\ntrain_ts = model_df.loc[\n    model_df.index <= train_end,\n    \"total_sales\"\n].copy()\n\nvalidation_ts = model_df.loc[\n    (model_df.index > train_end) &\n    (model_df.index <= validation_end),\n    \"total_sales\"\n].copy()\n\ntest_ts = model_df.loc[\n    model_df.index >= test_start,\n    \"total_sales\"\n].copy()\n\n\n# ------------------------------------------------------------\n# 4. Display split information\n# ------------------------------------------------------------\n\nprint(\"=\" * 70)\nprint(\"TIME SERIES SPLIT\")\nprint(\"=\" * 70)\n\nprint(\"\\nTRAIN\")\nprint(\"Start:\", train_ts.index.min())\nprint(\"End:\", train_ts.index.max())\nprint(\"Observations:\", len(train_ts))\n\nprint(\"\\nVALIDATION\")\nprint(\"Start:\", validation_ts.index.min())\nprint(\"End:\", validation_ts.index.max())\nprint(\"Observations:\", len(validation_ts))\n\nprint(\"\\nTEST\")\nprint(\"Start:\", test_ts.index.min())\nprint(\"End:\", test_ts.index.max())\nprint(\"Observations:\", len(test_ts))\n\n\n# ------------------------------------------------------------\n# 5. Check for missing values\n# ------------------------------------------------------------\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"MISSING VALUE CHECK\")\nprint(\"=\" * 70)\n\nprint(\"Train missing:\", train_ts.isna().sum())\nprint(\"Validation missing:\", validation_ts.isna().sum())\nprint(\"Test missing:\", test_ts.isna().sum())\n\n\n# ------------------------------------------------------------\n# 6. Check that dates do not overlap\n# ------------------------------------------------------------\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"OVERLAP CHECK\")\nprint(\"=\" * 70)\n\nprint(\n    \"Train ∩ Validation:\",\n    len(train_ts.index.intersection(validation_ts.index))\n)\n\nprint(\n    \"Validation ∩ Test:\",\n    len(validation_ts.index.intersection(test_ts.index))\n)\n\nprint(\n    \"Train ∩ Test:\",\n    len(train_ts.index.intersection(test_ts.index))\n)\n\n\n# ------------------------------------------------------------\n# 7. Plot the split\n# ------------------------------------------------------------\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    train_ts.index,\n    train_ts.values,\n    label=\"Train\"\n)\n\nplt.plot(\n    validation_ts.index,\n    validation_ts.values,\n    label=\"Validation\"\n)\n\nplt.plot(\n    test_ts.index,\n    test_ts.values,\n    label=\"Test\"\n)\n\nplt.axvline(\n    train_end,\n    linestyle=\"--\"\n)\n\nplt.axvline(\n    validation_end,\n    linestyle=\"--\"\n)\n\nplt.title(\"Chronological Train / Validation / Test Split\")\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\nplt.xticks(rotation=45)\n\nplt.tight_layout()\nplt.show()\n\n\n# ------------------------------------------------------------\n# 8. Save the splits\n# ------------------------------------------------------------\n\ntrain_ts.to_csv(\n    \"/kaggle/working/train_timeseries.csv\",\n    header=True\n)\n\nvalidation_ts.to_csv(\n    \"/kaggle/working/validation_timeseries.csv\",\n    header=True\n)\n\ntest_ts.to_csv(\n    \"/kaggle/working/test_timeseries.csv\",\n    header=True\n)\n\nprint(\"\\nSaved:\")\nprint(\"/kaggle/working/train_timeseries.csv\")\nprint(\"/kaggle/working/validation_timeseries.csv\")\nprint(\"/kaggle/working/test_timeseries.csv\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:21:42.378483Z","iopub.execute_input":"2026-09-17T10:21:42.378808Z","iopub.status.idle":"2026-09-17T10:21:42.748929Z","shell.execute_reply.started":"2026-09-17T10:21:42.378775Z","shell.execute_reply":"2026-09-17T10:21:42.747514Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# STEP 19 — BASELINE FORECASTS + EVALUATION METRICS\n# ============================================================\n\nfrom sklearn.metrics import (\n    mean_absolute_error,\n    mean_squared_error\n)\n\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n\n\n# ============================================================\n# 1. FUNCTION TO CALCULATE FORECAST METRICS\n# ============================================================\n\ndef calculate_metrics(actual, predicted, name=\"Model\"):\n\n    actual = pd.Series(actual).reset_index(drop=True)\n    predicted = pd.Series(predicted).reset_index(drop=True)\n\n    # -------------------------------\n    # MAE\n    # -------------------------------\n    mae = mean_absolute_error(\n        actual,\n        predicted\n    )\n\n    # -------------------------------\n    # RMSE\n    # -------------------------------\n    rmse = np.sqrt(\n        mean_squared_error(\n            actual,\n            predicted\n        )\n    )\n\n    # -------------------------------\n    # MAPE\n    # Exclude actual zero values\n    # -------------------------------\n    non_zero_mask = actual != 0\n\n    mape = (\n        np.mean(\n            np.abs(\n                (\n                    actual[non_zero_mask]\n                    - predicted[non_zero_mask]\n                )\n                /\n                actual[non_zero_mask]\n            )\n        )\n        * 100\n    )\n\n    excluded_from_mape = (\n        (~non_zero_mask).sum()\n    )\n\n    # -------------------------------\n    # Directional Accuracy\n    # -------------------------------\n    # Compare:\n    # actual direction:\n    # actual[t] - actual[t-1]\n    #\n    # predicted direction:\n    # predicted[t] - actual[t-1]\n    #\n    # This measures whether the forecast\n    # gets the direction of movement correct.\n    # -------------------------------\n\n    actual_direction = np.sign(\n        actual.diff()\n    )\n\n    predicted_direction = np.sign(\n        predicted - actual.shift(1)\n    )\n\n    direction_mask = (\n        actual_direction.notna()\n        & predicted_direction.notna()\n    )\n\n    directional_accuracy = (\n        np.mean(\n            actual_direction[direction_mask]\n            ==\n            predicted_direction[direction_mask]\n        )\n        * 100\n    )\n\n    return {\n        \"Model\": name,\n        \"MAE\": mae,\n        \"RMSE\": rmse,\n        \"MAPE (%)\": mape,\n        \"Directional Accuracy (%)\":\n            directional_accuracy,\n        \"MAPE Zero Exclusions\":\n            excluded_from_mape\n    }\n\n\n# ============================================================\n# 2. NAIVE FORECAST\n# ============================================================\n\n# Every validation prediction is the\n# last known training value.\n\nnaive_predictions = pd.Series(\n    train_ts.iloc[-1],\n    index=validation_ts.index\n)\n\n\n# ============================================================\n# 3. SEASONAL NAIVE FORECAST\n# ============================================================\n\n# Use the value from 7 days earlier.\n\ncombined_series = pd.concat([\n    train_ts,\n    validation_ts\n])\n\nseasonal_naive_predictions = (\n    combined_series\n    .shift(7)\n    .loc[validation_ts.index]\n)\n\n\n# ============================================================\n# 4. EVALUATE NAIVE\n# ============================================================\n\nnaive_metrics = calculate_metrics(\n    validation_ts,\n    naive_predictions,\n    \"Naive\"\n)\n\n\n# ============================================================\n# 5. EVALUATE SEASONAL NAIVE\n# ============================================================\n\nseasonal_naive_metrics = calculate_metrics(\n    validation_ts,\n    seasonal_naive_predictions,\n    \"Seasonal Naive (7-day)\"\n)\n\n\n# ============================================================\n# 6. CREATE MODEL COMPARISON TABLE\n# ============================================================\n\nbaseline_results = pd.DataFrame([\n    naive_metrics,\n    seasonal_naive_metrics\n])\n\nprint(\"=\" * 70)\nprint(\"BASELINE MODEL RESULTS — VALIDATION SET\")\nprint(\"=\" * 70)\n\ndisplay(\n    baseline_results.round(3)\n)\n\n\n# ============================================================\n# 7. COMPARE ACTUAL VS NAIVE\n# ============================================================\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    validation_ts.index,\n    validation_ts.values,\n    label=\"Actual\"\n)\n\nplt.plot(\n    validation_ts.index,\n    naive_predictions.values,\n    label=\"Naive Forecast\"\n)\n\nplt.title(\n    \"Actual vs Naive Forecast — Validation Set\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\nplt.xticks(rotation=45)\n\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 8. COMPARE ACTUAL VS SEASONAL NAIVE\n# ============================================================\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    validation_ts.index,\n    validation_ts.values,\n    label=\"Actual\"\n)\n\nplt.plot(\n    validation_ts.index,\n    seasonal_naive_predictions.values,\n    label=\"Seasonal Naive (7-day)\"\n)\n\nplt.title(\n    \"Actual vs Seasonal Naive Forecast — Validation Set\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\nplt.xticks(rotation=45)\n\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 9. SAVE BASELINE RESULTS\n# ============================================================\n\nbaseline_results.to_csv(\n    \"/kaggle/working/baseline_results.csv\",\n    index=False\n)\n\nprint(\"\\nSaved:\")\nprint(\"/kaggle/working/baseline_results.csv\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:22:10.035545Z","iopub.execute_input":"2026-09-17T10:22:10.035892Z","iopub.status.idle":"2026-09-17T10:22:11.220913Z","shell.execute_reply.started":"2026-09-17T10:22:10.035865Z","shell.execute_reply":"2026-09-17T10:22:11.219586Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# STEP 20 — ARIMA MODEL\n# ============================================================\n\nfrom statsmodels.tsa.arima.model import ARIMA\nimport warnings\n\nwarnings.filterwarnings(\"ignore\")\n\n\n# ============================================================\n# 1. TRAIN ARIMA(2,1,2)\n# ============================================================\n\nprint(\"=\" * 70)\nprint(\"TRAINING ARIMA(2,1,2)\")\nprint(\"=\" * 70)\n\narima_model = ARIMA(\n    train_ts,\n    order=(2, 1, 2)\n)\n\narima_result = arima_model.fit()\n\nprint(arima_result.summary())\n\n\n# ============================================================\n# 2. FORECAST VALIDATION PERIOD\n# ============================================================\n\narima_predictions = arima_result.forecast(\n    steps=len(validation_ts)\n)\n\n# Make sure prediction index matches validation dates\narima_predictions.index = validation_ts.index\n\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"ARIMA FORECAST\")\nprint(\"=\" * 70)\n\nprint(\"Forecast observations:\", len(arima_predictions))\n\ndisplay(\n    arima_predictions.head(10)\n)\n\n\n# ============================================================\n# 3. EVALUATE ARIMA\n# ============================================================\n\narima_metrics = calculate_metrics(\n    validation_ts,\n    arima_predictions,\n    \"ARIMA(2,1,2)\"\n)\n\narima_results = pd.DataFrame([\n    arima_metrics\n])\n\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"ARIMA RESULTS — VALIDATION SET\")\nprint(\"=\" * 70)\n\ndisplay(\n    arima_results.round(3)\n)\n\n\n# ============================================================\n# 4. COMBINE WITH BASELINES\n# ============================================================\n\nall_model_results = pd.concat(\n    [\n        baseline_results,\n        arima_results\n    ],\n    ignore_index=True\n)\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"BASELINE vs ARIMA\")\nprint(\"=\" * 70)\n\ndisplay(\n    all_model_results.round(3)\n)\n\n\n# ============================================================\n# 5. ACTUAL VS ARIMA FORECAST\n# ============================================================\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    validation_ts.index,\n    validation_ts.values,\n    label=\"Actual\"\n)\n\nplt.plot(\n    arima_predictions.index,\n    arima_predictions.values,\n    label=\"ARIMA(2,1,2)\"\n)\n\nplt.title(\n    \"Actual vs ARIMA(2,1,2) — Validation Set\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\nplt.xticks(rotation=45)\n\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 6. ACTUAL vs ARIMA vs SEASONAL NAIVE\n# ============================================================\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    validation_ts.index,\n    validation_ts.values,\n    label=\"Actual\"\n)\n\nplt.plot(\n    validation_ts.index,\n    seasonal_naive_predictions.values,\n    label=\"Seasonal Naive (7-day)\"\n)\n\nplt.plot(\n    validation_ts.index,\n    arima_predictions.values,\n    label=\"ARIMA(2,1,2)\"\n)\n\nplt.title(\n    \"Actual vs Seasonal Naive vs ARIMA\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\nplt.xticks(rotation=45)\n\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 7. SAVE ARIMA RESULTS\n# ============================================================\n\narima_results.to_csv(\n    \"/kaggle/working/arima_results.csv\",\n    index=False\n)\n\narima_predictions.to_csv(\n    \"/kaggle/working/arima_validation_predictions.csv\",\n    header=[\"prediction\"]\n)\n\nprint(\"\\nSaved:\")\nprint(\"/kaggle/working/arima_results.csv\")\nprint(\"/kaggle/working/arima_validation_predictions.csv\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:22:41.352777Z","iopub.execute_input":"2026-09-17T10:22:41.353269Z","iopub.status.idle":"2026-09-17T10:22:43.809787Z","shell.execute_reply.started":"2026-09-17T10:22:41.35322Z","shell.execute_reply":"2026-09-17T10:22:43.808476Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# STEP 21 — SARIMA WITH 7-DAY SEASONALITY\n# ============================================================\n\nfrom statsmodels.tsa.statespace.sarimax import SARIMAX\nimport warnings\n\nwarnings.filterwarnings(\"ignore\")\n\n\n# ============================================================\n# 1. TRAIN SARIMA\n# ============================================================\n\nprint(\"=\" * 70)\nprint(\"TRAINING SARIMA(1,1,1)(1,0,1,7)\")\nprint(\"=\" * 70)\n\nsarima_model = SARIMAX(\n    train_ts,\n    order=(1, 1, 1),\n    seasonal_order=(1, 0, 1, 7),\n    enforce_stationarity=False,\n    enforce_invertibility=False\n)\n\nsarima_result = sarima_model.fit(\n    disp=False\n)\n\nprint(sarima_result.summary())\n\n\n# ============================================================\n# 2. FORECAST VALIDATION PERIOD\n# ============================================================\n\nsarima_predictions = sarima_result.forecast(\n    steps=len(validation_ts)\n)\n\nsarima_predictions.index = validation_ts.index\n\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"SARIMA FORECAST\")\nprint(\"=\" * 70)\n\nprint(\"Forecast observations:\",\n      len(sarima_predictions))\n\ndisplay(\n    sarima_predictions.head(10)\n)\n\n\n# ============================================================\n# 3. EVALUATE SARIMA\n# ============================================================\n\nsarima_metrics = calculate_metrics(\n    validation_ts,\n    sarima_predictions,\n    \"SARIMA(1,1,1)(1,0,1,7)\"\n)\n\nsarima_results = pd.DataFrame([\n    sarima_metrics\n])\n\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"SARIMA RESULTS — VALIDATION SET\")\nprint(\"=\" * 70)\n\ndisplay(\n    sarima_results.round(3)\n)\n\n\n# ============================================================\n# 4. COMPARE ALL CURRENT MODELS\n# ============================================================\n\nall_model_results = pd.concat(\n    [\n        baseline_results,\n        arima_results,\n        sarima_results\n    ],\n    ignore_index=True\n)\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"CURRENT MODEL COMPARISON\")\nprint(\"=\" * 70)\n\ndisplay(\n    all_model_results.round(3)\n)\n\n\n# ============================================================\n# 5. ACTUAL VS SARIMA\n# ============================================================\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    validation_ts.index,\n    validation_ts.values,\n    label=\"Actual\"\n)\n\nplt.plot(\n    sarima_predictions.index,\n    sarima_predictions.values,\n    label=\"SARIMA(1,1,1)(1,0,1,7)\"\n)\n\nplt.title(\n    \"Actual vs SARIMA — Validation Set\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\nplt.xticks(rotation=45)\n\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 6. ACTUAL VS SEASONAL NAIVE VS SARIMA\n# ============================================================\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    validation_ts.index,\n    validation_ts.values,\n    label=\"Actual\"\n)\n\nplt.plot(\n    validation_ts.index,\n    seasonal_naive_predictions.values,\n    label=\"Seasonal Naive (7-day)\"\n)\n\nplt.plot(\n    validation_ts.index,\n    sarima_predictions.values,\n    label=\"SARIMA\"\n)\n\nplt.title(\n    \"Actual vs Seasonal Naive vs SARIMA\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\nplt.xticks(rotation=45)\n\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 7. SAVE RESULTS\n# ============================================================\n\nsarima_results.to_csv(\n    \"/kaggle/working/sarima_results.csv\",\n    index=False\n)\n\nsarima_predictions.to_csv(\n    \"/kaggle/working/sarima_validation_predictions.csv\",\n    header=[\"prediction\"]\n)\n\nprint(\"\\nSaved:\")\nprint(\"/kaggle/working/sarima_results.csv\")\nprint(\"/kaggle/working/sarima_validation_predictions.csv\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:23:16.516871Z","iopub.execute_input":"2026-09-17T10:23:16.517285Z","iopub.status.idle":"2026-09-17T10:23:18.595819Z","shell.execute_reply.started":"2026-09-17T10:23:16.517256Z","shell.execute_reply":"2026-09-17T10:23:18.594777Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# STEP 22 — EXPONENTIAL SMOOTHING / HOLT-WINTERS\n# ============================================================\n\nfrom statsmodels.tsa.holtwinters import ExponentialSmoothing\nimport warnings\n\nwarnings.filterwarnings(\"ignore\")\n\n\n# ============================================================\n# 1. TRAIN HOLT-WINTERS MODEL\n# ============================================================\n\nprint(\"=\" * 70)\nprint(\"TRAINING EXPONENTIAL SMOOTHING\")\nprint(\"=\" * 70)\n\nexp_model = ExponentialSmoothing(\n    train_ts,\n    trend=\"add\",\n    seasonal=\"add\",\n    seasonal_periods=7,\n    initialization_method=\"estimated\"\n)\n\nexp_result = exp_model.fit(\n    optimized=True\n)\n\nprint(exp_result.summary())\n\n\n# ============================================================\n# 2. FORECAST VALIDATION PERIOD\n# ============================================================\n\nexp_predictions = exp_result.forecast(\n    len(validation_ts)\n)\n\nexp_predictions.index = validation_ts.index\n\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"EXPONENTIAL SMOOTHING FORECAST\")\nprint(\"=\" * 70)\n\nprint(\"Forecast observations:\",\n      len(exp_predictions))\n\ndisplay(\n    exp_predictions.head(10)\n)\n\n\n# ============================================================\n# 3. EVALUATE MODEL\n# ============================================================\n\nexp_metrics = calculate_metrics(\n    validation_ts,\n    exp_predictions,\n    \"Exponential Smoothing\"\n)\n\nexp_results = pd.DataFrame([\n    exp_metrics\n])\n\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"EXPONENTIAL SMOOTHING RESULTS\")\nprint(\"=\" * 70)\n\ndisplay(\n    exp_results.round(3)\n)\n\n\n# ============================================================\n# 4. CURRENT MODEL COMPARISON\n# ============================================================\n\nall_model_results = pd.concat(\n    [\n        baseline_results,\n        arima_results,\n        sarima_results,\n        exp_results\n    ],\n    ignore_index=True\n)\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"CURRENT MODEL COMPARISON\")\nprint(\"=\" * 70)\n\ndisplay(\n    all_model_results.round(3)\n)\n\n\n# ============================================================\n# 5. ACTUAL VS EXPONENTIAL SMOOTHING\n# ============================================================\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    validation_ts.index,\n    validation_ts.values,\n    label=\"Actual\"\n)\n\nplt.plot(\n    exp_predictions.index,\n    exp_predictions.values,\n    label=\"Exponential Smoothing\"\n)\n\nplt.title(\n    \"Actual vs Exponential Smoothing — Validation Set\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\nplt.xticks(rotation=45)\n\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 6. ACTUAL VS SARIMA VS EXPONENTIAL SMOOTHING\n# ============================================================\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    validation_ts.index,\n    validation_ts.values,\n    label=\"Actual\"\n)\n\nplt.plot(\n    sarima_predictions.index,\n    sarima_predictions.values,\n    label=\"SARIMA\"\n)\n\nplt.plot(\n    exp_predictions.index,\n    exp_predictions.values,\n    label=\"Exponential Smoothing\"\n)\n\nplt.title(\n    \"Actual vs SARIMA vs Exponential Smoothing\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\nplt.xticks(rotation=45)\n\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 7. SAVE RESULTS\n# ============================================================\n\nexp_results.to_csv(\n    \"/kaggle/working/exponential_smoothing_results.csv\",\n    index=False\n)\n\nexp_predictions.to_csv(\n    \"/kaggle/working/exponential_smoothing_predictions.csv\",\n    header=[\"prediction\"]\n)\n\nprint(\"\\nSaved:\")\nprint(\n    \"/kaggle/working/exponential_smoothing_results.csv\"\n)\n\nprint(\n    \"/kaggle/working/exponential_smoothing_predictions.csv\"\n)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:23:52.73891Z","iopub.execute_input":"2026-09-17T10:23:52.739296Z","iopub.status.idle":"2026-09-17T10:23:53.877456Z","shell.execute_reply.started":"2026-09-17T10:23:52.739265Z","shell.execute_reply":"2026-09-17T10:23:53.876297Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# STEP 23 — PROPHET FORECASTING MODEL\n# ============================================================\n\nfrom prophet import Prophet\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport warnings\n\nwarnings.filterwarnings(\"ignore\")\n\n\n# ============================================================\n# 1. PREPARE DATA FOR PROPHET\n# ============================================================\n\nprophet_train = (\n    train_ts\n    .reset_index()\n    .rename(\n        columns={\n            \"date\": \"ds\",\n            \"total_sales\": \"y\"\n        }\n    )\n)\n\nprint(\"=\" * 70)\nprint(\"PROPHET TRAINING DATA\")\nprint(\"=\" * 70)\n\nprint(\"Shape:\", prophet_train.shape)\n\ndisplay(prophet_train.head())\n\n\n# ============================================================\n# 2. CREATE PROPHET MODEL\n# ============================================================\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"TRAINING PROPHET\")\nprint(\"=\" * 70)\n\nprophet_model = Prophet(\n    yearly_seasonality=True,\n    weekly_seasonality=True,\n    daily_seasonality=False,\n    seasonality_mode=\"additive\",\n    changepoint_prior_scale=0.05\n)\n\n\n# ============================================================\n# 3. TRAIN MODEL\n# ============================================================\n\nprophet_model.fit(prophet_train)\n\n\n# ============================================================\n# 4. CREATE FUTURE DATES\n# ============================================================\n\nfuture = prophet_model.make_future_dataframe(\n    periods=len(validation_ts),\n    freq=\"D\",\n    include_history=False\n)\n\nprint(\"\\nFuture dates:\")\ndisplay(future.head())\ndisplay(future.tail())\n\n\n# ============================================================\n# 5. GENERATE FORECAST\n# ============================================================\n\nprophet_forecast = prophet_model.predict(\n    future\n)\n\n\n# ============================================================\n# 6. EXTRACT PREDICTIONS\n# ============================================================\n\nprophet_predictions = (\n    prophet_forecast[\"yhat\"]\n    .copy()\n)\n\nprophet_predictions.index = validation_ts.index\n\n\n# ============================================================\n# 7. EVALUATE PROPHET\n# ============================================================\n\nprophet_metrics = calculate_metrics(\n    validation_ts,\n    prophet_predictions,\n    \"Prophet\"\n)\n\nprophet_results = pd.DataFrame([\n    prophet_metrics\n])\n\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"PROPHET RESULTS — VALIDATION SET\")\nprint(\"=\" * 70)\n\ndisplay(\n    prophet_results.round(3)\n)\n\n\n# ============================================================\n# 8. CURRENT MODEL COMPARISON\n# ============================================================\n\nall_model_results = pd.concat(\n    [\n        baseline_results,\n        arima_results,\n        sarima_results,\n        exp_results,\n        prophet_results\n    ],\n    ignore_index=True\n)\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"CURRENT MODEL COMPARISON\")\nprint(\"=\" * 70)\n\ndisplay(\n    all_model_results.round(3)\n)\n\n\n# ============================================================\n# 9. ACTUAL VS PROPHET\n# ============================================================\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    validation_ts.index,\n    validation_ts.values,\n    label=\"Actual\"\n)\n\nplt.plot(\n    prophet_predictions.index,\n    prophet_predictions.values,\n    label=\"Prophet\"\n)\n\nplt.title(\n    \"Actual vs Prophet — Validation Set\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\nplt.xticks(rotation=45)\n\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 10. ACTUAL VS ALL MAIN MODELS\n# ============================================================\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    validation_ts.index,\n    validation_ts.values,\n    label=\"Actual\"\n)\n\nplt.plot(\n    sarima_predictions.index,\n    sarima_predictions.values,\n    label=\"SARIMA\"\n)\n\nplt.plot(\n    exp_predictions.index,\n    exp_predictions.values,\n    label=\"Exponential Smoothing\"\n)\n\nplt.plot(\n    prophet_predictions.index,\n    prophet_predictions.values,\n    label=\"Prophet\"\n)\n\nplt.title(\n    \"Actual vs SARIMA vs Exponential Smoothing vs Prophet\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\nplt.xticks(rotation=45)\n\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 11. PROPHET FORECAST COMPONENTS\n# ============================================================\n\nfig1 = prophet_model.plot_components(\n    prophet_model.predict(future)\n)\n\nplt.show()\n\n\n# ============================================================\n# 12. SAVE RESULTS\n# ============================================================\n\nprophet_results.to_csv(\n    \"/kaggle/working/prophet_results.csv\",\n    index=False\n)\n\nprophet_predictions.to_csv(\n    \"/kaggle/working/prophet_validation_predictions.csv\",\n    header=[\"prediction\"]\n)\n\nprint(\"\\nSaved:\")\nprint(\"/kaggle/working/prophet_results.csv\")\nprint(\"/kaggle/working/prophet_validation_predictions.csv\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:24:18.069531Z","iopub.execute_input":"2026-09-17T10:24:18.070589Z","iopub.status.idle":"2026-09-17T10:24:22.457679Z","shell.execute_reply.started":"2026-09-17T10:24:18.070528Z","shell.execute_reply":"2026-09-17T10:24:22.456635Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# STEP 24 — LSTM FORECASTING MODEL\n# ============================================================\n\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport warnings\n\nfrom sklearn.preprocessing import MinMaxScaler\nfrom tensorflow.keras.models import Sequential\nfrom tensorflow.keras.layers import LSTM, Dense, Dropout\nfrom tensorflow.keras.callbacks import EarlyStopping\n\nwarnings.filterwarnings(\"ignore\")\n\n\n# ============================================================\n# 1. PREPARE TRAINING DATA\n# ============================================================\n\ntrain_values = train_ts.values.reshape(-1, 1)\n\nprint(\"=\" * 70)\nprint(\"LSTM DATA PREPARATION\")\nprint(\"=\" * 70)\n\nprint(\"Training observations:\", len(train_values))\n\n\n# ============================================================\n# 2. SCALE TRAINING DATA\n# ============================================================\n\nscaler = MinMaxScaler(\n    feature_range=(0, 1)\n)\n\ntrain_scaled = scaler.fit_transform(\n    train_values\n)\n\n\n# ============================================================\n# 3. CREATE SEQUENCES\n# ============================================================\n\nlookback = 30\n\nX_train = []\ny_train = []\n\nfor i in range(\n    lookback,\n    len(train_scaled)\n):\n    \n    X_train.append(\n        train_scaled[\n            i - lookback:i,\n            0\n        ]\n    )\n    \n    y_train.append(\n        train_scaled[i, 0]\n    )\n\nX_train = np.array(X_train)\ny_train = np.array(y_train)\n\n\n# LSTM expects:\n# samples × time_steps × features\n\nX_train = X_train.reshape(\n    X_train.shape[0],\n    X_train.shape[1],\n    1\n)\n\nprint(\"Lookback:\", lookback)\nprint(\"X_train shape:\", X_train.shape)\nprint(\"y_train shape:\", y_train.shape)\n\n\n# ============================================================\n# 4. BUILD LSTM MODEL\n# ============================================================\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"BUILDING LSTM MODEL\")\nprint(\"=\" * 70)\n\nlstm_model = Sequential([\n    \n    LSTM(\n        64,\n        return_sequences=True,\n        input_shape=(lookback, 1)\n    ),\n    \n    Dropout(0.2),\n    \n    LSTM(\n        32,\n        return_sequences=False\n    ),\n    \n    Dropout(0.2),\n    \n    Dense(16, activation=\"relu\"),\n    \n    Dense(1)\n])\n\n\nlstm_model.compile(\n    optimizer=\"adam\",\n    loss=\"mse\"\n)\n\nlstm_model.summary()\n\n\n# ============================================================\n# 5. TRAIN LSTM\n# ============================================================\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"TRAINING LSTM\")\nprint(\"=\" * 70)\n\nearly_stopping = EarlyStopping(\n    monitor=\"val_loss\",\n    patience=8,\n    restore_best_weights=True\n)\n\nhistory = lstm_model.fit(\n    X_train,\n    y_train,\n    epochs=50,\n    batch_size=32,\n    validation_split=0.1,\n    shuffle=False,\n    callbacks=[early_stopping],\n    verbose=1\n)\n\n\n# ============================================================\n# 6. PLOT TRAINING HISTORY\n# ============================================================\n\nplt.figure(figsize=(10, 5))\n\nplt.plot(\n    history.history[\"loss\"],\n    label=\"Training Loss\"\n)\n\nplt.plot(\n    history.history[\"val_loss\"],\n    label=\"Validation Loss\"\n)\n\nplt.title(\"LSTM Training History\")\nplt.xlabel(\"Epoch\")\nplt.ylabel(\"MSE Loss\")\n\nplt.legend()\n\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 7. RECURSIVE FORECASTING\n# ============================================================\n\n# Start with the final 30 training observations\n# and recursively predict the validation period.\n\nlast_sequence = train_scaled[\n    -lookback:\n].flatten().tolist()\n\nlstm_predictions_scaled = []\n\nfor _ in range(len(validation_ts)):\n    \n    input_sequence = np.array(\n        last_sequence[-lookback:]\n    )\n    \n    input_sequence = input_sequence.reshape(\n        1,\n        lookback,\n        1\n    )\n    \n    next_prediction = (\n        lstm_model\n        .predict(\n            input_sequence,\n            verbose=0\n        )[0, 0]\n    )\n    \n    lstm_predictions_scaled.append(\n        next_prediction\n    )\n    \n    # Add prediction to sequence\n    last_sequence.append(\n        next_prediction\n    )\n\n\n# ============================================================\n# 8. CONVERT BACK TO ORIGINAL SCALE\n# ============================================================\n\nlstm_predictions = scaler.inverse_transform(\n    np.array(\n        lstm_predictions_scaled\n    ).reshape(-1, 1)\n).flatten()\n\nlstm_predictions = pd.Series(\n    lstm_predictions,\n    index=validation_ts.index\n)\n\n\n# ============================================================\n# 9. EVALUATE LSTM\n# ============================================================\n\nlstm_metrics = calculate_metrics(\n    validation_ts,\n    lstm_predictions,\n    \"LSTM\"\n)\n\nlstm_results = pd.DataFrame([\n    lstm_metrics\n])\n\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"LSTM RESULTS — VALIDATION SET\")\nprint(\"=\" * 70)\n\ndisplay(\n    lstm_results.round(3)\n)\n\n\n# ============================================================\n# 10. FINAL MODEL COMPARISON\n# ============================================================\n\nall_model_results = pd.concat(\n    [\n        baseline_results,\n        arima_results,\n        sarima_results,\n        exp_results,\n        prophet_results,\n        lstm_results\n    ],\n    ignore_index=True\n)\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"ALL MODEL COMPARISON\")\nprint(\"=\" * 70)\n\ndisplay(\n    all_model_results.round(3)\n)\n\n\n# ============================================================\n# 11. SORT BY MAE\n# ============================================================\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"MODELS SORTED BY MAE\")\nprint(\"=\" * 70)\n\ndisplay(\n    all_model_results\n    .sort_values(\"MAE\")\n    .reset_index(drop=True)\n    .round(3)\n)\n\n\n# ============================================================\n# 12. ACTUAL VS LSTM\n# ============================================================\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    validation_ts.index,\n    validation_ts.values,\n    label=\"Actual\"\n)\n\nplt.plot(\n    lstm_predictions.index,\n    lstm_predictions.values,\n    label=\"LSTM\"\n)\n\nplt.title(\n    \"Actual vs LSTM — Validation Set\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\nplt.xticks(rotation=45)\n\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 13. ACTUAL VS ALL MAJOR MODELS\n# ============================================================\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    validation_ts.index,\n    validation_ts.values,\n    label=\"Actual\"\n)\n\nplt.plot(\n    sarima_predictions.index,\n    sarima_predictions.values,\n    label=\"SARIMA\"\n)\n\nplt.plot(\n    exp_predictions.index,\n    exp_predictions.values,\n    label=\"Exponential Smoothing\"\n)\n\nplt.plot(\n    prophet_predictions.index,\n    prophet_predictions.values,\n    label=\"Prophet\"\n)\n\nplt.plot(\n    lstm_predictions.index,\n    lstm_predictions.values,\n    label=\"LSTM\"\n)\n\nplt.title(\n    \"Actual vs Major Forecasting Models\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\nplt.xticks(rotation=45)\n\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 14. SAVE RESULTS\n# ============================================================\n\nlstm_results.to_csv(\n    \"/kaggle/working/lstm_results.csv\",\n    index=False\n)\n\nlstm_predictions.to_csv(\n    \"/kaggle/working/lstm_validation_predictions.csv\",\n    header=[\"prediction\"]\n)\n\nlstm_model.save(\n    \"/kaggle/working/lstm_demand_model.keras\"\n)\n\nprint(\"\\nSaved:\")\nprint(\"/kaggle/working/lstm_results.csv\")\nprint(\"/kaggle/working/lstm_validation_predictions.csv\")\nprint(\"/kaggle/working/lstm_demand_model.keras\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:24:56.561098Z","iopub.execute_input":"2026-09-17T10:24:56.561496Z","iopub.status.idle":"2026-09-17T10:25:53.589802Z","shell.execute_reply.started":"2026-09-17T10:24:56.561467Z","shell.execute_reply":"2026-09-17T10:25:53.58829Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# STEP 25 — SARIMA HYPERPARAMETER TUNING\n# ============================================================\n\nimport itertools\nimport warnings\nimport numpy as np\nimport pandas as pd\nfrom statsmodels.tsa.statespace.sarimax import SARIMAX\n\nwarnings.filterwarnings(\"ignore\")\n\n\n# ============================================================\n# 1. DEFINE SARIMA PARAMETER GRID\n# ============================================================\n\n# Non-seasonal parameters\np_values = [0, 1, 2]\nd_values = [1]\nq_values = [0, 1, 2]\n\n# Seasonal parameters\nP_values = [0, 1]\nD_values = [0]\nQ_values = [0, 1]\n\nseasonal_period = 7\n\n\n# ============================================================\n# 2. GENERATE PARAMETER COMBINATIONS\n# ============================================================\n\nparameter_combinations = list(\n    itertools.product(\n        p_values,\n        d_values,\n        q_values,\n        P_values,\n        D_values,\n        Q_values\n    )\n)\n\nprint(\"=\" * 70)\nprint(\"SARIMA HYPERPARAMETER TUNING\")\nprint(\"=\" * 70)\n\nprint(\n    \"Total combinations:\",\n    len(parameter_combinations)\n)\n\n\n# ============================================================\n# 3. GRID SEARCH\n# ============================================================\n\ntuning_results = []\n\nfor i, params in enumerate(parameter_combinations, start=1):\n\n    p, d, q, P, D, Q = params\n\n    order = (p, d, q)\n\n    seasonal_order = (\n        P,\n        D,\n        Q,\n        seasonal_period\n    )\n\n    print(\n        f\"[{i}/{len(parameter_combinations)}] \"\n        f\"SARIMA{order}x{seasonal_order}\"\n    )\n\n    try:\n\n        model = SARIMAX(\n            train_ts,\n            order=order,\n            seasonal_order=seasonal_order,\n            enforce_stationarity=False,\n            enforce_invertibility=False\n        )\n\n        result = model.fit(\n            disp=False\n        )\n\n        predictions = result.forecast(\n            steps=len(validation_ts)\n        )\n\n        predictions.index = validation_ts.index\n\n        metrics = calculate_metrics(\n            validation_ts,\n            predictions,\n            f\"SARIMA{order}x{seasonal_order}\"\n        )\n\n        metrics[\"AIC\"] = result.aic\n        metrics[\"Order\"] = str(order)\n        metrics[\"Seasonal Order\"] = str(\n            seasonal_order\n        )\n\n        tuning_results.append(metrics)\n\n    except Exception as e:\n\n        print(\n            \"   Failed:\",\n            str(e)[:100]\n        )\n\n\n# ============================================================\n# 4. CREATE RESULTS TABLE\n# ============================================================\n\ntuning_results_df = pd.DataFrame(\n    tuning_results\n)\n\n\n# ============================================================\n# 5. SORT BY MAE\n# ============================================================\n\ntuning_results_by_mae = (\n    tuning_results_df\n    .sort_values(\"MAE\")\n    .reset_index(drop=True)\n)\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"SARIMA MODELS SORTED BY MAE\")\nprint(\"=\" * 70)\n\ndisplay(\n    tuning_results_by_mae[\n        [\n            \"Model\",\n            \"MAE\",\n            \"RMSE\",\n            \"MAPE (%)\",\n            \"Directional Accuracy (%)\",\n            \"AIC\",\n            \"Order\",\n            \"Seasonal Order\"\n        ]\n    ].round(3)\n)\n\n\n# ============================================================\n# 6. SORT BY RMSE\n# ============================================================\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"SARIMA MODELS SORTED BY RMSE\")\nprint(\"=\" * 70)\n\ndisplay(\n    tuning_results_df\n    .sort_values(\"RMSE\")\n    .reset_index(drop=True)[\n        [\n            \"Model\",\n            \"MAE\",\n            \"RMSE\",\n            \"MAPE (%)\",\n            \"Directional Accuracy (%)\",\n            \"AIC\",\n            \"Order\",\n            \"Seasonal Order\"\n        ]\n    ].round(3)\n)\n\n\n# ============================================================\n# 7. SORT BY MAPE\n# ============================================================\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"SARIMA MODELS SORTED BY MAPE\")\nprint(\"=\" * 70)\n\ndisplay(\n    tuning_results_df\n    .sort_values(\"MAPE (%)\")\n    .reset_index(drop=True)[\n        [\n            \"Model\",\n            \"MAE\",\n            \"RMSE\",\n            \"MAPE (%)\",\n            \"Directional Accuracy (%)\",\n            \"AIC\",\n            \"Order\",\n            \"Seasonal Order\"\n        ]\n    ].round(3)\n)\n\n\n# ============================================================\n# 8. SELECT BEST MODEL BY MAE\n# ============================================================\n\nbest_sarima_row = (\n    tuning_results_by_mae.iloc[0]\n)\n\nbest_order = eval(\n    best_sarima_row[\"Order\"]\n)\n\nbest_seasonal_order = eval(\n    best_sarima_row[\"Seasonal Order\"]\n)\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"BEST SARIMA MODEL\")\nprint(\"=\" * 70)\n\nprint(\n    \"Order:\",\n    best_order\n)\n\nprint(\n    \"Seasonal Order:\",\n    best_seasonal_order\n)\n\nprint(\n    \"MAE:\",\n    round(best_sarima_row[\"MAE\"], 3)\n)\n\nprint(\n    \"RMSE:\",\n    round(best_sarima_row[\"RMSE\"], 3)\n)\n\nprint(\n    \"MAPE:\",\n    round(best_sarima_row[\"MAPE (%)\"], 3)\n)\n\n\n# ============================================================\n# 9. TRAIN THE BEST MODEL AGAIN\n# ============================================================\n\nbest_sarima_model = SARIMAX(\n    train_ts,\n    order=best_order,\n    seasonal_order=best_seasonal_order,\n    enforce_stationarity=False,\n    enforce_invertibility=False\n)\n\nbest_sarima_result = (\n    best_sarima_model.fit(\n        disp=False\n    )\n)\n\n\n# ============================================================\n# 10. GENERATE BEST MODEL VALIDATION FORECAST\n# ============================================================\n\nbest_sarima_predictions = (\n    best_sarima_result\n    .forecast(\n        steps=len(validation_ts)\n    )\n)\n\nbest_sarima_predictions.index = (\n    validation_ts.index\n)\n\n\n# ============================================================\n# 11. FINAL TUNED MODEL METRICS\n# ============================================================\n\nbest_sarima_metrics = calculate_metrics(\n    validation_ts,\n    best_sarima_predictions,\n    \"Tuned SARIMA\"\n)\n\nbest_sarima_metrics_df = pd.DataFrame([\n    best_sarima_metrics\n])\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"TUNED SARIMA VALIDATION RESULTS\")\nprint(\"=\" * 70)\n\ndisplay(\n    best_sarima_metrics_df.round(3)\n)\n\n\n# ============================================================\n# 12. ORIGINAL vs TUNED SARIMA\n# ============================================================\n\noriginal_sarima_mae = (\n    sarima_results[\"MAE\"].iloc[0]\n)\n\ntuned_sarima_mae = (\n    best_sarima_metrics_df[\"MAE\"].iloc[0]\n)\n\nimprovement = (\n    (original_sarima_mae - tuned_sarima_mae)\n    / original_sarima_mae\n) * 100\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"TUNING IMPROVEMENT\")\nprint(\"=\" * 70)\n\nprint(\n    \"Original SARIMA MAE:\",\n    round(original_sarima_mae, 3)\n)\n\nprint(\n    \"Tuned SARIMA MAE:\",\n    round(tuned_sarima_mae, 3)\n)\n\nprint(\n    \"MAE improvement (%):\",\n    round(improvement, 2)\n)\n\n\n# ============================================================\n# 13. PLOT TUNED SARIMA\n# ============================================================\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    validation_ts.index,\n    validation_ts.values,\n    label=\"Actual\"\n)\n\nplt.plot(\n    best_sarima_predictions.index,\n    best_sarima_predictions.values,\n    label=\"Tuned SARIMA\"\n)\n\nplt.title(\n    \"Actual vs Tuned SARIMA — Validation Set\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\nplt.xticks(rotation=45)\n\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 14. SAVE TUNING RESULTS\n# ============================================================\n\ntuning_results_df.to_csv(\n    \"/kaggle/working/sarima_tuning_results.csv\",\n    index=False\n)\n\nbest_sarima_predictions.to_csv(\n    \"/kaggle/working/tuned_sarima_validation_predictions.csv\",\n    header=[\"prediction\"]\n)\n\nprint(\"\\nSaved:\")\nprint(\"/kaggle/working/sarima_tuning_results.csv\")\nprint(\n    \"/kaggle/working/\"\n    \"tuned_sarima_validation_predictions.csv\"\n)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:26:46.390689Z","iopub.execute_input":"2026-09-17T10:26:46.392223Z","iopub.status.idle":"2026-09-17T10:27:22.618151Z","shell.execute_reply.started":"2026-09-17T10:26:46.39218Z","shell.execute_reply":"2026-09-17T10:27:22.616585Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# STEP 26 — EXPANDING-WINDOW ROLLING BACKTESTING\n# ============================================================\n\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom statsmodels.tsa.statespace.sarimax import SARIMAX\nimport warnings\n\nwarnings.filterwarnings(\"ignore\")\n\n\n# ============================================================\n# 1. BACKTEST SETTINGS\n# ============================================================\n\nBACKTEST_HORIZON = 90\nMIN_TRAIN_SIZE = 730\nSTEP_SIZE = 90\n\n# Tuned model obtained in Step 25\nBACKTEST_ORDER = (2, 1, 2)\nBACKTEST_SEASONAL_ORDER = (1, 0, 1, 7)\n\n\n# ============================================================\n# 2. STORE RESULTS\n# ============================================================\n\nbacktest_results = []\nbacktest_predictions = []\n\n\n# ============================================================\n# 3. CREATE EXPANDING-WINDOW FOLDS\n# ============================================================\n\nfold_number = 1\ntrain_size = MIN_TRAIN_SIZE\n\nwhile train_size + BACKTEST_HORIZON <= len(train_ts):\n\n    fold_train = train_ts.iloc[:train_size]\n\n    fold_actual = train_ts.iloc[\n        train_size:\n        train_size + BACKTEST_HORIZON\n    ]\n\n    print(\"=\" * 70)\n    print(f\"BACKTEST FOLD {fold_number}\")\n    print(\"=\" * 70)\n\n    print(\"Training:\")\n    print(\n        fold_train.index.min(),\n        \"to\",\n        fold_train.index.max()\n    )\n\n    print(\"Forecast:\")\n    print(\n        fold_actual.index.min(),\n        \"to\",\n        fold_actual.index.max()\n    )\n\n    print(\"Training observations:\", len(fold_train))\n    print(\"Forecast observations:\", len(fold_actual))\n\n\n    # --------------------------------------------------------\n    # Train SARIMA\n    # --------------------------------------------------------\n\n    sarima_model = SARIMAX(\n        fold_train,\n        order=BACKTEST_ORDER,\n        seasonal_order=BACKTEST_SEASONAL_ORDER,\n        enforce_stationarity=False,\n        enforce_invertibility=False\n    )\n\n    sarima_result = sarima_model.fit(\n        disp=False\n    )\n\n\n    # --------------------------------------------------------\n    # Forecast\n    # --------------------------------------------------------\n\n    sarima_pred = sarima_result.forecast(\n        steps=BACKTEST_HORIZON\n    )\n\n    sarima_pred.index = fold_actual.index\n\n\n    # --------------------------------------------------------\n    # Seasonal Naive baseline\n    # --------------------------------------------------------\n\n    combined_fold = pd.concat([\n        fold_train,\n        fold_actual\n    ])\n\n    seasonal_naive_pred = (\n        combined_fold\n        .shift(7)\n        .loc[fold_actual.index]\n    )\n\n\n    # --------------------------------------------------------\n    # Calculate metrics\n    # --------------------------------------------------------\n\n    sarima_metrics = calculate_metrics(\n        fold_actual,\n        sarima_pred,\n        f\"SARIMA Fold {fold_number}\"\n    )\n\n    seasonal_naive_metrics = calculate_metrics(\n        fold_actual,\n        seasonal_naive_pred,\n        f\"Seasonal Naive Fold {fold_number}\"\n    )\n\n    sarima_metrics[\"Fold\"] = fold_number\n    seasonal_naive_metrics[\"Fold\"] = fold_number\n\n    sarima_metrics[\"Forecast Start\"] = (\n        fold_actual.index.min()\n    )\n\n    seasonal_naive_metrics[\"Forecast Start\"] = (\n        fold_actual.index.min()\n    )\n\n\n    backtest_results.append(\n        sarima_metrics\n    )\n\n    backtest_results.append(\n        seasonal_naive_metrics\n    )\n\n\n    # --------------------------------------------------------\n    # Save predictions\n    # --------------------------------------------------------\n\n    fold_predictions = pd.DataFrame({\n        \"date\": fold_actual.index,\n        \"actual\": fold_actual.values,\n        \"sarima_prediction\": sarima_pred.values,\n        \"seasonal_naive_prediction\":\n            seasonal_naive_pred.values,\n        \"fold\": fold_number\n    })\n\n    backtest_predictions.append(\n        fold_predictions\n    )\n\n\n    fold_number += 1\n    train_size += STEP_SIZE\n\n\n# ============================================================\n# 4. COMBINE BACKTEST RESULTS\n# ============================================================\n\nbacktest_results_df = pd.DataFrame(\n    backtest_results\n)\n\nbacktest_predictions_df = pd.concat(\n    backtest_predictions,\n    ignore_index=True\n)\n\n\n# ============================================================\n# 5. DISPLAY FOLD RESULTS\n# ============================================================\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"BACKTEST RESULTS BY FOLD\")\nprint(\"=\" * 70)\n\ndisplay(\n    backtest_results_df[\n        [\n            \"Fold\",\n            \"Model\",\n            \"MAE\",\n            \"RMSE\",\n            \"MAPE (%)\",\n            \"Directional Accuracy (%)\"\n        ]\n    ].round(3)\n)\n\n\n# ============================================================\n# 6. AVERAGE BACKTEST PERFORMANCE\n# ============================================================\n\naverage_backtest = (\n    backtest_results_df\n    .groupby(\"Model\")[\n        [\n            \"MAE\",\n            \"RMSE\",\n            \"MAPE (%)\",\n            \"Directional Accuracy (%)\"\n        ]\n    ]\n    .mean()\n    .reset_index()\n)\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"AVERAGE BACKTEST PERFORMANCE\")\nprint(\"=\" * 70)\n\ndisplay(\n    average_backtest.round(3)\n)\n\n\n# ============================================================\n# 7. MEDIAN BACKTEST PERFORMANCE\n# ============================================================\n\nmedian_backtest = (\n    backtest_results_df\n    .groupby(\"Model\")[\n        [\n            \"MAE\",\n            \"RMSE\",\n            \"MAPE (%)\",\n            \"Directional Accuracy (%)\"\n        ]\n    ]\n    .median()\n    .reset_index()\n)\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"MEDIAN BACKTEST PERFORMANCE\")\nprint(\"=\" * 70)\n\ndisplay(\n    median_backtest.round(3)\n)\n\n\n# ============================================================\n# 8. PLOT BACKTEST FOLD ERRORS\n# ============================================================\n\nsarima_fold_results = (\n    backtest_results_df[\n        backtest_results_df[\"Model\"].str.startswith(\"SARIMA\")\n    ]\n    .copy()\n)\n\nseasonal_fold_results = (\n    backtest_results_df[\n        backtest_results_df[\"Model\"].str.startswith(\"Seasonal Naive\")\n    ]\n    .copy()\n)\n\n\nplt.figure(figsize=(10, 5))\n\nplt.plot(\n    sarima_fold_results[\"Fold\"],\n    sarima_fold_results[\"MAE\"],\n    marker=\"o\",\n    label=\"Tuned SARIMA\"\n)\n\nplt.plot(\n    seasonal_fold_results[\"Fold\"],\n    seasonal_fold_results[\"MAE\"],\n    marker=\"o\",\n    label=\"Seasonal Naive\"\n)\n\nplt.title(\"Backtesting MAE by Fold\")\nplt.xlabel(\"Backtest Fold\")\nplt.ylabel(\"MAE\")\n\nplt.legend()\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 9. PLOT ONE BACKTEST FOLD\n# ============================================================\n\nfirst_fold = (\n    backtest_predictions_df[\n        backtest_predictions_df[\"fold\"] == 1\n    ]\n)\n\nplt.figure(figsize=(15, 6))\n\nplt.plot(\n    first_fold[\"date\"],\n    first_fold[\"actual\"],\n    label=\"Actual\"\n)\n\nplt.plot(\n    first_fold[\"date\"],\n    first_fold[\"sarima_prediction\"],\n    label=\"Tuned SARIMA\"\n)\n\nplt.plot(\n    first_fold[\"date\"],\n    first_fold[\"seasonal_naive_prediction\"],\n    label=\"Seasonal Naive\"\n)\n\nplt.title(\"Backtest Fold 1 — Actual vs Forecast\")\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\nplt.xticks(rotation=45)\n\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 10. SAVE BACKTEST RESULTS\n# ============================================================\n\nbacktest_results_df.to_csv(\n    \"/kaggle/working/backtest_results.csv\",\n    index=False\n)\n\naverage_backtest.to_csv(\n    \"/kaggle/working/backtest_average_results.csv\",\n    index=False\n)\n\nbacktest_predictions_df.to_csv(\n    \"/kaggle/working/backtest_predictions.csv\",\n    index=False\n)\n\nprint(\"\\nSaved:\")\nprint(\"/kaggle/working/backtest_results.csv\")\nprint(\"/kaggle/working/backtest_average_results.csv\")\nprint(\"/kaggle/working/backtest_predictions.csv\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:28:35.451636Z","iopub.execute_input":"2026-09-17T10:28:35.452593Z","iopub.status.idle":"2026-09-17T10:28:51.504751Z","shell.execute_reply.started":"2026-09-17T10:28:35.452555Z","shell.execute_reply":"2026-09-17T10:28:51.503724Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# STEP 26A — CORRECT BACKTEST SUMMARY\n# ============================================================\n\n# ------------------------------------------------------------\n# 1. Create a clean model category\n# ------------------------------------------------------------\n\nbacktest_results_df[\"Model Type\"] = (\n    backtest_results_df[\"Model\"]\n    .str.replace(\n        r\" Fold \\d+\",\n        \"\",\n        regex=True\n    )\n)\n\n\n# ============================================================\n# 2. AVERAGE PERFORMANCE ACROSS FOLDS\n# ============================================================\n\ncorrect_average_backtest = (\n    backtest_results_df\n    .groupby(\"Model Type\")[\n        [\n            \"MAE\",\n            \"RMSE\",\n            \"MAPE (%)\",\n            \"Directional Accuracy (%)\"\n        ]\n    ]\n    .agg([\"mean\", \"std\"])\n    .reset_index()\n)\n\nprint(\"=\" * 70)\nprint(\"CORRECT AVERAGE BACKTEST PERFORMANCE\")\nprint(\"=\" * 70)\n\ndisplay(\n    correct_average_backtest.round(3)\n)\n\n\n# ============================================================\n# 3. MEDIAN PERFORMANCE ACROSS FOLDS\n# ============================================================\n\ncorrect_median_backtest = (\n    backtest_results_df\n    .groupby(\"Model Type\")[\n        [\n            \"MAE\",\n            \"RMSE\",\n            \"MAPE (%)\",\n            \"Directional Accuracy (%)\"\n        ]\n    ]\n    .median()\n    .reset_index()\n)\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"CORRECT MEDIAN BACKTEST PERFORMANCE\")\nprint(\"=\" * 70)\n\ndisplay(\n    correct_median_backtest.round(3)\n)\n\n\n# ============================================================\n# 4. PERFORMANCE BY FOLD\n# ============================================================\n\nfold_comparison = (\n    backtest_results_df[\n        [\n            \"Fold\",\n            \"Model Type\",\n            \"MAE\",\n            \"RMSE\",\n            \"MAPE (%)\",\n            \"Directional Accuracy (%)\"\n        ]\n    ]\n    .sort_values(\n        [\"Fold\", \"Model Type\"]\n    )\n)\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"PERFORMANCE BY FOLD\")\nprint(\"=\" * 70)\n\ndisplay(\n    fold_comparison.round(3)\n)\n\n\n# ============================================================\n# 5. COUNT HOW MANY FOLDS SARIMA BEATS SEASONAL NAIVE\n# ============================================================\n\nsarima_fold = (\n    backtest_results_df[\n        backtest_results_df[\"Model Type\"]\n        == \"SARIMA\"\n    ][\n        [\n            \"Fold\",\n            \"MAE\",\n            \"RMSE\",\n            \"MAPE (%)\",\n            \"Directional Accuracy (%)\"\n        ]\n    ]\n    .set_index(\"Fold\")\n)\n\nseasonal_fold = (\n    backtest_results_df[\n        backtest_results_df[\"Model Type\"]\n        == \"Seasonal Naive\"\n    ][\n        [\n            \"Fold\",\n            \"MAE\",\n            \"RMSE\",\n            \"MAPE (%)\",\n            \"Directional Accuracy (%)\"\n        ]\n    ]\n    .set_index(\"Fold\")\n)\n\n\nsarima_beats_mae = (\n    sarima_fold[\"MAE\"]\n    < seasonal_fold[\"MAE\"]\n)\n\nsarima_beats_rmse = (\n    sarima_fold[\"RMSE\"]\n    < seasonal_fold[\"RMSE\"]\n)\n\nsarima_beats_mape = (\n    sarima_fold[\"MAPE (%)\"]\n    < seasonal_fold[\"MAPE (%)\"]\n)\n\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"SARIMA VS SEASONAL NAIVE — FOLD WIN COUNT\")\nprint(\"=\" * 70)\n\nprint(\n    \"SARIMA better MAE:\",\n    sarima_beats_mae.sum(),\n    \"out of\",\n    len(sarima_beats_mae)\n)\n\nprint(\n    \"SARIMA better RMSE:\",\n    sarima_beats_rmse.sum(),\n    \"out of\",\n    len(sarima_beats_rmse)\n)\n\nprint(\n    \"SARIMA better MAPE:\",\n    sarima_beats_mape.sum(),\n    \"out of\",\n    len(sarima_beats_mape)\n)\n\n\n# ============================================================\n# 6. POOLED BACKTEST METRICS\n# ============================================================\n\npooled_mae = {}\npooled_rmse = {}\npooled_mape = {}\npooled_directional = {}\n\nfor model_type in [\"SARIMA\", \"Seasonal Naive\"]:\n\n    subset = backtest_predictions_df.copy()\n\n    if model_type == \"SARIMA\":\n\n        actual = subset[\"actual\"]\n        predicted = subset[\"sarima_prediction\"]\n\n    else:\n\n        actual = subset[\"actual\"]\n        predicted = subset[\n            \"seasonal_naive_prediction\"\n        ]\n\n\n    actual = pd.Series(actual).reset_index(drop=True)\n    predicted = pd.Series(predicted).reset_index(drop=True)\n\n\n    pooled_mae[model_type] = mean_absolute_error(\n        actual,\n        predicted\n    )\n\n    pooled_rmse[model_type] = np.sqrt(\n        mean_squared_error(\n            actual,\n            predicted\n        )\n    )\n\n\n    non_zero = actual != 0\n\n    pooled_mape[model_type] = (\n        np.mean(\n            np.abs(\n                (\n                    actual[non_zero]\n                    - predicted[non_zero]\n                )\n                / actual[non_zero]\n            )\n        )\n        * 100\n    )\n\n\n    actual_direction = np.sign(\n        actual.diff()\n    )\n\n    predicted_direction = np.sign(\n        predicted - actual.shift(1)\n    )\n\n    valid_direction = (\n        actual_direction.notna()\n        &\n        predicted_direction.notna()\n    )\n\n    pooled_directional[model_type] = (\n        np.mean(\n            actual_direction[valid_direction]\n            ==\n            predicted_direction[valid_direction]\n        )\n        * 100\n    )\n\n\npooled_results = pd.DataFrame({\n    \"Model\": [\n        \"SARIMA\",\n        \"Seasonal Naive\"\n    ],\n    \"MAE\": [\n        pooled_mae[\"SARIMA\"],\n        pooled_mae[\"Seasonal Naive\"]\n    ],\n    \"RMSE\": [\n        pooled_rmse[\"SARIMA\"],\n        pooled_rmse[\"Seasonal Naive\"]\n    ],\n    \"MAPE (%)\": [\n        pooled_mape[\"SARIMA\"],\n        pooled_mape[\"Seasonal Naive\"]\n    ],\n    \"Directional Accuracy (%)\": [\n        pooled_directional[\"SARIMA\"],\n        pooled_directional[\"Seasonal Naive\"]\n    ]\n})\n\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"POOLED BACKTEST PERFORMANCE\")\nprint(\"=\" * 70)\n\ndisplay(\n    pooled_results.round(3)\n)\n\n\n# ============================================================\n# 7. SAVE CORRECTED RESULTS\n# ============================================================\n\ncorrect_average_backtest.to_csv(\n    \"/kaggle/working/\"\n    \"correct_average_backtest_results.csv\",\n    index=False\n)\n\ncorrect_median_backtest.to_csv(\n    \"/kaggle/working/\"\n    \"correct_median_backtest_results.csv\",\n    index=False\n)\n\npooled_results.to_csv(\n    \"/kaggle/working/\"\n    \"pooled_backtest_results.csv\",\n    index=False\n)\n\nprint(\"\\nSaved corrected backtest summaries.\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:29:53.237297Z","iopub.execute_input":"2026-09-17T10:29:53.238394Z","iopub.status.idle":"2026-09-17T10:29:53.343829Z","shell.execute_reply.started":"2026-09-17T10:29:53.238356Z","shell.execute_reply":"2026-09-17T10:29:53.342424Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# STEP 27 — FINAL TEST-SET EVALUATION\n# ============================================================\n\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n\nfrom statsmodels.tsa.statespace.sarimax import SARIMAX\n\n\n# ============================================================\n# 1. TRAIN TUNED SARIMA ON TRAIN + VALIDATION\n# ============================================================\n\n# We have already used validation for model tuning.\n# Now combine train + validation to give the selected\n# SARIMA model all historical information before the\n# untouched test period.\n\npre_test_series = pd.concat([\n    train_ts,\n    validation_ts\n])\n\nprint(\"=\" * 70)\nprint(\"FINAL TEST EVALUATION\")\nprint(\"=\" * 70)\n\nprint(\"Training period:\")\nprint(\n    pre_test_series.index.min(),\n    \"to\",\n    pre_test_series.index.max()\n)\n\nprint(\"\\nTest period:\")\nprint(\n    test_ts.index.min(),\n    \"to\",\n    test_ts.index.max()\n)\n\nprint(\"\\nTraining observations:\",\n      len(pre_test_series))\n\nprint(\"Test observations:\",\n      len(test_ts))\n\n\n# ============================================================\n# 2. FIT TUNED SARIMA\n# ============================================================\n\nfinal_test_model = SARIMAX(\n    pre_test_series,\n    order=(2, 1, 2),\n    seasonal_order=(1, 0, 1, 7),\n    enforce_stationarity=False,\n    enforce_invertibility=False\n)\n\nfinal_test_result = final_test_model.fit(\n    disp=False\n)\n\n\n# ============================================================\n# 3. FORECAST TEST PERIOD\n# ============================================================\n\ntest_predictions = final_test_result.forecast(\n    steps=len(test_ts)\n)\n\ntest_predictions.index = test_ts.index\n\n\n# ============================================================\n# 4. SEASONAL NAIVE TEST FORECAST\n# ============================================================\n\ncombined_pretest = pd.concat([\n    pre_test_series,\n    test_ts\n])\n\nseasonal_naive_test_predictions = (\n    combined_pretest\n    .shift(7)\n    .loc[test_ts.index]\n)\n\n\n# ============================================================\n# 5. EVALUATE TUNED SARIMA\n# ============================================================\n\nfinal_sarima_test_metrics = calculate_metrics(\n    test_ts,\n    test_predictions,\n    \"Tuned SARIMA\"\n)\n\n\n# ============================================================\n# 6. EVALUATE SEASONAL NAIVE\n# ============================================================\n\nfinal_seasonal_naive_test_metrics = calculate_metrics(\n    test_ts,\n    seasonal_naive_test_predictions,\n    \"Seasonal Naive\"\n)\n\n\n# ============================================================\n# 7. TEST RESULTS TABLE\n# ============================================================\n\nfinal_test_results = pd.DataFrame([\n    final_sarima_test_metrics,\n    final_seasonal_naive_test_metrics\n])\n\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"FINAL TEST RESULTS\")\nprint(\"=\" * 70)\n\ndisplay(\n    final_test_results.round(3)\n)\n\n\n# ============================================================\n# 8. COMPARE WITH VALIDATION RESULTS\n# ============================================================\n\ncomparison = pd.DataFrame({\n    \"Model\": [\n        \"Tuned SARIMA\",\n        \"Seasonal Naive\"\n    ],\n\n    \"Validation MAE\": [\n        best_sarima_metrics_df[\"MAE\"].iloc[0],\n        baseline_results[\n            baseline_results[\"Model\"]\n            == \"Seasonal Naive (7-day)\"\n        ][\"MAE\"].iloc[0]\n    ],\n\n    \"Test MAE\": [\n        final_sarima_test_metrics[\"MAE\"],\n        final_seasonal_naive_test_metrics[\"MAE\"]\n    ],\n\n    \"Validation RMSE\": [\n        best_sarima_metrics_df[\"RMSE\"].iloc[0],\n        baseline_results[\n            baseline_results[\"Model\"]\n            == \"Seasonal Naive (7-day)\"\n        ][\"RMSE\"].iloc[0]\n    ],\n\n    \"Test RMSE\": [\n        final_sarima_test_metrics[\"RMSE\"],\n        final_seasonal_naive_test_metrics[\"RMSE\"]\n    ],\n\n    \"Validation MAPE (%)\": [\n        best_sarima_metrics_df[\"MAPE (%)\"].iloc[0],\n        baseline_results[\n            baseline_results[\"Model\"]\n            == \"Seasonal Naive (7-day)\"\n        ][\"MAPE (%)\"].iloc[0]\n    ],\n\n    \"Test MAPE (%)\": [\n        final_sarima_test_metrics[\"MAPE (%)\"],\n        final_seasonal_naive_test_metrics[\"MAPE (%)\"]\n    ],\n\n    \"Validation Direction (%)\": [\n        best_sarima_metrics_df[\n            \"Directional Accuracy (%)\"\n        ].iloc[0],\n\n        baseline_results[\n            baseline_results[\"Model\"]\n            == \"Seasonal Naive (7-day)\"\n        ][\"Directional Accuracy (%)\"].iloc[0]\n    ],\n\n    \"Test Direction (%)\": [\n        final_sarima_test_metrics[\n            \"Directional Accuracy (%)\"\n        ],\n\n        final_seasonal_naive_test_metrics[\n            \"Directional Accuracy (%)\"\n        ]\n    ]\n})\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"VALIDATION VS TEST\")\nprint(\"=\" * 70)\n\ndisplay(\n    comparison.round(3)\n)\n\n\n# ============================================================\n# 9. PLOT TEST FORECAST\n# ============================================================\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    test_ts.index,\n    test_ts.values,\n    label=\"Actual\"\n)\n\nplt.plot(\n    test_predictions.index,\n    test_predictions.values,\n    label=\"Tuned SARIMA\"\n)\n\nplt.plot(\n    seasonal_naive_test_predictions.index,\n    seasonal_naive_test_predictions.values,\n    label=\"Seasonal Naive\"\n)\n\nplt.title(\n    \"Final Test Set — Actual vs Forecasts\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\nplt.xticks(rotation=45)\n\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 10. SAVE FINAL TEST RESULTS\n# ============================================================\n\nfinal_test_results.to_csv(\n    \"/kaggle/working/final_test_results.csv\",\n    index=False\n)\n\ncomparison.to_csv(\n    \"/kaggle/working/validation_vs_test_comparison.csv\",\n    index=False\n)\n\ntest_predictions.to_csv(\n    \"/kaggle/working/final_sarima_test_predictions.csv\",\n    header=[\"prediction\"]\n)\n\nprint(\"\\nSaved:\")\nprint(\"/kaggle/working/final_test_results.csv\")\nprint(\"/kaggle/working/validation_vs_test_comparison.csv\")\nprint(\"/kaggle/working/final_sarima_test_predictions.csv\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:30:24.712372Z","iopub.execute_input":"2026-09-17T10:30:24.713575Z","iopub.status.idle":"2026-09-17T10:30:33.217637Z","shell.execute_reply.started":"2026-09-17T10:30:24.713535Z","shell.execute_reply":"2026-09-17T10:30:33.216543Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# STEP 28 — FINAL 6-MONTH AND 12-MONTH FORECAST\n# ============================================================\n\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n\nfrom statsmodels.tsa.statespace.sarimax import SARIMAX\n\n\n# ============================================================\n# 1. FINAL HISTORICAL SERIES\n# ============================================================\n\n# Use ALL available historical data\n# from 2013-01-01 through 2017-08-15.\n\nfinal_series = (\n    forecast_df\n    .set_index(\"date\")[\"total_sales\"]\n    .sort_index()\n    .asfreq(\"D\")\n)\n\nprint(\"=\" * 70)\nprint(\"FINAL FORECAST MODEL\")\nprint(\"=\" * 70)\n\nprint(\"Historical start:\",\n      final_series.index.min())\n\nprint(\"Historical end:\",\n      final_series.index.max())\n\nprint(\"Historical observations:\",\n      len(final_series))\n\n\n# ============================================================\n# 2. TRAIN FINAL TUNED SARIMA\n# ============================================================\n\nfinal_model = SARIMAX(\n    final_series,\n\n    order=(2, 1, 2),\n\n    seasonal_order=(1, 0, 1, 7),\n\n    enforce_stationarity=False,\n\n    enforce_invertibility=False\n)\n\nfinal_result = final_model.fit(\n    disp=False\n)\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"FINAL SARIMA MODEL\")\nprint(\"=\" * 70)\n\nprint(\"Order:\",\n      (2, 1, 2))\n\nprint(\"Seasonal Order:\",\n      (1, 0, 1, 7))\n\nprint(\"\\nModel fitted successfully.\")\n\n\n# ============================================================\n# 3. CREATE 12-MONTH FORECAST\n# ============================================================\n\nFORECAST_DAYS = 365\n\nfinal_forecast_result = (\n    final_result.get_forecast(\n        steps=FORECAST_DAYS\n    )\n)\n\n\n# ============================================================\n# 4. EXTRACT POINT FORECAST\n# ============================================================\n\nforecast_mean = (\n    final_forecast_result.predicted_mean\n)\n\nforecast_mean.name = \"forecast\"\n\n\n# ============================================================\n# 5. EXTRACT 95% CONFIDENCE INTERVAL\n# ============================================================\n\nconfidence_intervals = (\n    final_forecast_result\n    .conf_int(alpha=0.05)\n)\n\nconfidence_intervals.columns = [\n    \"lower_95\",\n    \"upper_95\"\n]\n\n\n# ============================================================\n# 6. CREATE FINAL FORECAST DATAFRAME\n# ============================================================\n\nforecast_output = pd.concat(\n    [\n        forecast_mean,\n        confidence_intervals\n    ],\n    axis=1\n)\n\nforecast_output.index.name = \"date\"\n\n# Demand cannot be negative from a business perspective.\n# Keep the statistical SARIMA forecast unchanged,\n# but create an inventory-use forecast clipped at zero.\n\nforecast_output[\"inventory_forecast\"] = (\n    forecast_output[\"forecast\"]\n    .clip(lower=0)\n)\n\nforecast_output[\"lower_95_inventory\"] = (\n    forecast_output[\"lower_95\"]\n    .clip(lower=0)\n)\n\nforecast_output[\"upper_95_inventory\"] = (\n    forecast_output[\"upper_95\"]\n    .clip(lower=0)\n)\n\n\n# ============================================================\n# 7. ADD FORECAST HORIZON\n# ============================================================\n\nforecast_output[\"days_ahead\"] = (\n    np.arange(1, FORECAST_DAYS + 1)\n)\n\nforecast_output[\"forecast_month\"] = (\n    forecast_output.index.to_period(\"M\")\n    .astype(str)\n)\n\n\n# ============================================================\n# 8. DISPLAY FIRST 10 FORECAST DAYS\n# ============================================================\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"FIRST 10 FORECAST DAYS\")\nprint(\"=\" * 70)\n\ndisplay(\n    forecast_output.head(10).round(2)\n)\n\n\n# ============================================================\n# 9. DISPLAY LAST 10 FORECAST DAYS\n# ============================================================\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"LAST 10 FORECAST DAYS\")\nprint(\"=\" * 70)\n\ndisplay(\n    forecast_output.tail(10).round(2)\n)\n\n\n# ============================================================\n# 10. 6-MONTH FORECAST\n# ============================================================\n\nforecast_6m = forecast_output.iloc[:182].copy()\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"6-MONTH FORECAST\")\nprint(\"=\" * 70)\n\nprint(\"Start:\",\n      forecast_6m.index.min())\n\nprint(\"End:\",\n      forecast_6m.index.max())\n\nprint(\"Average daily forecast:\",\n      round(\n          forecast_6m[\"inventory_forecast\"].mean(),\n          2\n      ))\n\nprint(\"Total forecast demand:\",\n      round(\n          forecast_6m[\"inventory_forecast\"].sum(),\n          2\n      ))\n\n\n# ============================================================\n# 11. 12-MONTH FORECAST\n# ============================================================\n\nforecast_12m = forecast_output.copy()\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"12-MONTH FORECAST\")\nprint(\"=\" * 70)\n\nprint(\"Start:\",\n      forecast_12m.index.min())\n\nprint(\"End:\",\n      forecast_12m.index.max())\n\nprint(\"Average daily forecast:\",\n      round(\n          forecast_12m[\"inventory_forecast\"].mean(),\n          2\n      ))\n\nprint(\"Total forecast demand:\",\n      round(\n          forecast_12m[\"inventory_forecast\"].sum(),\n          2\n      ))\n\n\n# ============================================================\n# 12. MONTHLY FORECAST SUMMARY\n# ============================================================\n\nmonthly_forecast = (\n    forecast_output\n    .resample(\"ME\")\n    .agg(\n        forecast_units=(\n            \"inventory_forecast\",\n            \"sum\"\n        ),\n        average_daily_forecast=(\n            \"inventory_forecast\",\n            \"mean\"\n        ),\n        lower_95_units=(\n            \"lower_95_inventory\",\n            \"sum\"\n        ),\n        upper_95_units=(\n            \"upper_95_inventory\",\n            \"sum\"\n        )\n    )\n    .reset_index()\n)\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"MONTHLY FORECAST\")\nprint(\"=\" * 70)\n\ndisplay(\n    monthly_forecast.round(2)\n)\n\n\n# ============================================================\n# 13. FULL HISTORICAL + FORECAST VISUALIZATION\n# ============================================================\n\nplt.figure(figsize=(18, 7))\n\n# Historical\nplt.plot(\n    final_series.index,\n    final_series.values,\n    label=\"Historical Sales\"\n)\n\n# Forecast\nplt.plot(\n    forecast_output.index,\n    forecast_output[\"inventory_forecast\"],\n    label=\"Forecast\"\n)\n\n# Confidence interval\nplt.fill_between(\n    forecast_output.index,\n    forecast_output[\"lower_95_inventory\"],\n    forecast_output[\"upper_95_inventory\"],\n    alpha=0.2,\n    label=\"95% Confidence Interval\"\n)\n\n# Forecast boundary\nplt.axvline(\n    final_series.index.max(),\n    linestyle=\"--\",\n    label=\"Forecast Start\"\n)\n\nplt.title(\n    \"Historical Sales and 12-Month SARIMA Forecast\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\nplt.xticks(rotation=45)\n\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 14. 6-MONTH FORECAST VISUALIZATION\n# ============================================================\n\nplt.figure(figsize=(16, 6))\n\n# Last 180 historical days\nrecent_history = final_series.tail(180)\n\nplt.plot(\n    recent_history.index,\n    recent_history.values,\n    label=\"Historical Sales\"\n)\n\nplt.plot(\n    forecast_6m.index,\n    forecast_6m[\"inventory_forecast\"],\n    label=\"6-Month Forecast\"\n)\n\nplt.fill_between(\n    forecast_6m.index,\n    forecast_6m[\"lower_95_inventory\"],\n    forecast_6m[\"upper_95_inventory\"],\n    alpha=0.2,\n    label=\"95% Confidence Interval\"\n)\n\nplt.axvline(\n    final_series.index.max(),\n    linestyle=\"--\",\n    label=\"Forecast Start\"\n)\n\nplt.title(\n    \"6-Month Demand Forecast with 95% Confidence Interval\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\nplt.xticks(rotation=45)\n\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 15. FORECAST UNCERTAINTY\n# ============================================================\n\nforecast_output[\"interval_width\"] = (\n    forecast_output[\"upper_95_inventory\"]\n    -\n    forecast_output[\"lower_95_inventory\"]\n)\n\nplt.figure(figsize=(15, 5))\n\nplt.plot(\n    forecast_output.index,\n    forecast_output[\"interval_width\"]\n)\n\nplt.title(\n    \"Forecast Uncertainty Over 12 Months\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"95% Confidence Interval Width\")\n\nplt.xticks(rotation=45)\n\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 16. SAVE FORECAST RESULTS\n# ============================================================\n\nforecast_output.to_csv(\n    \"/kaggle/working/final_12_month_forecast.csv\"\n)\n\nforecast_6m.to_csv(\n    \"/kaggle/working/final_6_month_forecast.csv\"\n)\n\nmonthly_forecast.to_csv(\n    \"/kaggle/working/monthly_forecast_summary.csv\",\n    index=False\n)\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"FILES SAVED\")\nprint(\"=\" * 70)\n\nprint(\n    \"/kaggle/working/final_12_month_forecast.csv\"\n)\n\nprint(\n    \"/kaggle/working/final_6_month_forecast.csv\"\n)\n\nprint(\n    \"/kaggle/working/monthly_forecast_summary.csv\"\n)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:31:14.94115Z","iopub.execute_input":"2026-09-17T10:31:14.941912Z","iopub.status.idle":"2026-09-17T10:31:20.153496Z","shell.execute_reply.started":"2026-09-17T10:31:14.941843Z","shell.execute_reply":"2026-09-17T10:31:20.152306Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# STEP 29 — INVENTORY / SUPPLY CHAIN RECOMMENDATIONS\n# ============================================================\n\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n\n\n# ============================================================\n# 1. CREATE DAILY INVENTORY PLANNING TABLE\n# ============================================================\n\ninventory_df = forecast_output.copy()\n\ninventory_df[\"forecast\"] = (\n    inventory_df[\"forecast\"].clip(lower=0)\n)\n\ninventory_df[\"lower_95\"] = (\n    inventory_df[\"lower_95\"].clip(lower=0)\n)\n\ninventory_df[\"upper_95\"] = (\n    inventory_df[\"upper_95\"].clip(lower=0)\n)\n\n\n# ------------------------------------------------------------\n# Forecast uncertainty\n# ------------------------------------------------------------\n\ninventory_df[\"uncertainty_range\"] = (\n    inventory_df[\"upper_95\"]\n    -\n    inventory_df[\"lower_95\"]\n)\n\n\ninventory_df[\"uncertainty_pct\"] = (\n    inventory_df[\"uncertainty_range\"]\n    /\n    inventory_df[\"forecast\"].replace(0, np.nan)\n) * 100\n\n\n# ============================================================\n# 2. SELECT PLANNING HORIZONS\n# ============================================================\n\nplanning_horizons = {\n    \"7_days\": 7,\n    \"14_days\": 14,\n    \"30_days\": 30\n}\n\n\nplanning_results = []\n\n\n# ============================================================\n# 3. CALCULATE INVENTORY PLANNING LEVELS\n# ============================================================\n\nfor horizon_name, days in planning_horizons.items():\n\n    horizon_data = inventory_df.iloc[:days]\n\n    expected_demand = (\n        horizon_data[\"forecast\"].sum()\n    )\n\n    conservative_demand = (\n        horizon_data[\"upper_95\"].sum()\n    )\n\n    lower_demand = (\n        horizon_data[\"lower_95\"].sum()\n    )\n\n    safety_buffer = (\n        conservative_demand\n        -\n        expected_demand\n    )\n\n    planning_results.append({\n\n        \"Planning Horizon\":\n            horizon_name,\n\n        \"Days\":\n            days,\n\n        \"Expected Demand\":\n            expected_demand,\n\n        \"Lower Demand Bound\":\n            lower_demand,\n\n        \"Upper Demand Bound\":\n            conservative_demand,\n\n        \"Safety Buffer\":\n            safety_buffer,\n\n        \"Recommended Planning Stock\":\n            conservative_demand\n    })\n\n\ninventory_planning = pd.DataFrame(\n    planning_results\n)\n\n\n# ============================================================\n# 4. MONTHLY FORECAST ANALYSIS\n# ============================================================\n\nmonthly_inventory = monthly_forecast.copy()\n\nmonthly_inventory[\"forecast_share\"] = (\n    monthly_inventory[\"forecast_units\"]\n    /\n    monthly_inventory[\"forecast_units\"].sum()\n    * 100\n)\n\nmonthly_inventory[\"uncertainty_width\"] = (\n    monthly_inventory[\"upper_95_units\"]\n    -\n    monthly_inventory[\"lower_95_units\"]\n)\n\n\n# ============================================================\n# 5. IDENTIFY HIGH AND LOW FORECAST MONTHS\n# ============================================================\n\ntop_months = (\n    monthly_inventory\n    .nlargest(\n        3,\n        \"forecast_units\"\n    )\n)\n\nlow_months = (\n    monthly_inventory\n    .nsmallest(\n        3,\n        \"forecast_units\"\n    )\n)\n\n\n# ============================================================\n# 6. PRINT INVENTORY PLANNING RESULTS\n# ============================================================\n\nprint(\"=\" * 70)\nprint(\"INVENTORY PLANNING — AGGREGATE DEMAND\")\nprint(\"=\" * 70)\n\ndisplay(\n    inventory_planning.round(2)\n)\n\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"TOP 3 HIGH-DEMAND FORECAST MONTHS\")\nprint(\"=\" * 70)\n\ndisplay(\n    top_months[\n        [\n            \"date\",\n            \"forecast_units\",\n            \"average_daily_forecast\",\n            \"forecast_share\"\n        ]\n    ].round(2)\n)\n\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"BOTTOM 3 FORECAST MONTHS\")\nprint(\"=\" * 70)\n\ndisplay(\n    low_months[\n        [\n            \"date\",\n            \"forecast_units\",\n            \"average_daily_forecast\",\n            \"forecast_share\"\n        ]\n    ].round(2)\n)\n\n\n# ============================================================\n# 7. 6-MONTH vs 12-MONTH PLANNING SUMMARY\n# ============================================================\n\nsix_month_summary = {\n    \"Period\": \"6 Months\",\n    \"Total Expected Demand\":\n        forecast_6m[\"forecast\"].sum(),\n    \"Average Daily Demand\":\n        forecast_6m[\"forecast\"].mean(),\n    \"Minimum Daily Forecast\":\n        forecast_6m[\"forecast\"].min(),\n    \"Maximum Daily Forecast\":\n        forecast_6m[\"forecast\"].max()\n}\n\ntwelve_month_summary = {\n    \"Period\": \"12 Months\",\n    \"Total Expected Demand\":\n        forecast_12m[\"forecast\"].sum(),\n    \"Average Daily Demand\":\n        forecast_12m[\"forecast\"].mean(),\n    \"Minimum Daily Forecast\":\n        forecast_12m[\"forecast\"].min(),\n    \"Maximum Daily Forecast\":\n        forecast_12m[\"forecast\"].max()\n}\n\nplanning_summary = pd.DataFrame([\n    six_month_summary,\n    twelve_month_summary\n])\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"6-MONTH VS 12-MONTH DEMAND SUMMARY\")\nprint(\"=\" * 70)\n\ndisplay(\n    planning_summary.round(2)\n)\n\n\n# ============================================================\n# 8. FORECAST UNCERTAINTY VISUALIZATION\n# ============================================================\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    inventory_df.index,\n    inventory_df[\"forecast\"],\n    label=\"Expected Demand\"\n)\n\nplt.fill_between(\n    inventory_df.index,\n    inventory_df[\"lower_95\"],\n    inventory_df[\"upper_95\"],\n    alpha=0.2,\n    label=\"95% Forecast Interval\"\n)\n\nplt.title(\n    \"Forecast Demand and Uncertainty for Inventory Planning\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\n\nplt.xticks(rotation=45)\n\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 9. MONTHLY DEMAND VISUALIZATION\n# ============================================================\n\nplt.figure(figsize=(14, 6))\n\nplt.bar(\n    monthly_inventory[\"date\"].astype(str),\n    monthly_inventory[\"forecast_units\"]\n)\n\nplt.title(\n    \"Forecast Monthly Demand\"\n)\n\nplt.xlabel(\"Month\")\nplt.ylabel(\"Forecast Unit Sales\")\n\nplt.xticks(\n    rotation=60\n)\n\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 10. GENERATE BUSINESS RECOMMENDATIONS\n# ============================================================\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"BUSINESS RECOMMENDATIONS\")\nprint(\"=\" * 70)\n\nprint(\"\"\"\n1. Inventory Planning\n   Use the forecast as the baseline demand plan and maintain\n   additional buffer stock during periods with wider uncertainty.\n\n2. Short-Term Replenishment\n   The 7-day and 14-day planning horizons can support frequent\n   replenishment decisions where supplier lead times are short.\n\n3. Longer Lead Times\n   For suppliers with approximately 30-day lead times, use the\n   30-day planning estimate together with the forecast uncertainty\n   buffer rather than relying only on point forecasts.\n\n4. Seasonal Planning\n   Increase procurement and warehouse capacity ahead of months\n   with higher forecast demand.\n\n5. Uncertainty Management\n   The confidence interval widens farther into the future.\n   Therefore, long-range procurement decisions should be reviewed\n   periodically and updated with newly observed sales.\n\n6. Product-Level Planning\n   These results represent aggregate demand. For operational SKU\n   inventory decisions, the same forecasting workflow should next\n   be applied at product/store level.\n\n7. Safety Stock\n   Treat the upper forecast bound as a conservative planning\n   scenario rather than guaranteed demand.\n\n8. Reforecasting\n   Refresh the forecast regularly as new sales observations become\n   available rather than treating the 12-month forecast as fixed.\n\"\"\")\n\n\n# ============================================================\n# 11. SAVE INVENTORY PLANNING OUTPUTS\n# ============================================================\n\ninventory_planning.to_csv(\n    \"/kaggle/working/inventory_planning.csv\",\n    index=False\n)\n\nmonthly_inventory.to_csv(\n    \"/kaggle/working/monthly_inventory_forecast.csv\",\n    index=False\n)\n\nplanning_summary.to_csv(\n    \"/kaggle/working/planning_summary.csv\",\n    index=False\n)\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"FILES SAVED\")\nprint(\"=\" * 70)\n\nprint(\"/kaggle/working/inventory_planning.csv\")\nprint(\"/kaggle/working/monthly_inventory_forecast.csv\")\nprint(\"/kaggle/working/planning_summary.csv\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:32:04.533081Z","iopub.execute_input":"2026-09-17T10:32:04.533469Z","iopub.status.idle":"2026-09-17T10:32:05.159111Z","shell.execute_reply.started":"2026-09-17T10:32:04.53344Z","shell.execute_reply":"2026-09-17T10:32:05.158099Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# STEP 30 — PRODUCT-LEVEL FORECASTING AND ACCURACY\n# ============================================================\n\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport warnings\n\nfrom statsmodels.tsa.statespace.sarimax import SARIMAX\n\nwarnings.filterwarnings(\"ignore\")\n\n\n# ============================================================\n# 1. FIND TOP 5 PRODUCTS BY HISTORICAL SALES\n# ============================================================\n\nprint(\"=\" * 70)\nprint(\"STEP 30 — FINDING TOP 5 PRODUCTS\")\nprint(\"=\" * 70)\n\nproduct_sales_parts = []\n\ntrain_reader = pd.read_csv(\n    train_path,\n    usecols=[\"item_nbr\", \"unit_sales\"],\n    chunksize=250000\n)\n\nfor chunk in train_reader:\n\n    product_totals = (\n        chunk\n        .groupby(\"item_nbr\")[\"unit_sales\"]\n        .sum()\n    )\n\n    product_sales_parts.append(\n        product_totals\n    )\n\n\nproduct_sales = (\n    pd.concat(product_sales_parts)\n    .groupby(level=0)\n    .sum()\n    .sort_values(ascending=False)\n)\n\n\ntop_5_products = (\n    product_sales\n    .head(5)\n    .reset_index()\n)\n\ntop_5_products.columns = [\n    \"item_nbr\",\n    \"historical_sales\"\n]\n\n\n# Add product information\ntop_5_products = top_5_products.merge(\n    items[\n        [\n            \"item_nbr\",\n            \"family\",\n            \"class\",\n            \"perishable\"\n        ]\n    ],\n    on=\"item_nbr\",\n    how=\"left\"\n)\n\n\nprint(\"\\nTop 5 products:\")\ndisplay(\n    top_5_products.round(2)\n)\n\n\n# ============================================================\n# 2. EXTRACT DAILY SALES FOR TOP 5 PRODUCTS\n# ============================================================\n\ntop_product_ids = (\n    top_5_products[\"item_nbr\"]\n    .tolist()\n)\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"EXTRACTING DAILY SALES FOR TOP 5 PRODUCTS\")\nprint(\"=\" * 70)\n\nproduct_daily_parts = []\n\ntrain_reader = pd.read_csv(\n    train_path,\n    usecols=[\n        \"date\",\n        \"item_nbr\",\n        \"unit_sales\"\n    ],\n    chunksize=250000\n)\n\nfor chunk in train_reader:\n\n    # Keep only top products\n    filtered = chunk[\n        chunk[\"item_nbr\"].isin(top_product_ids)\n    ].copy()\n\n    if len(filtered) == 0:\n        continue\n\n    filtered[\"date\"] = pd.to_datetime(\n        filtered[\"date\"]\n    )\n\n    grouped = (\n        filtered\n        .groupby(\n            [\"date\", \"item_nbr\"]\n        )[\"unit_sales\"]\n        .sum()\n        .reset_index()\n    )\n\n    product_daily_parts.append(\n        grouped\n    )\n\n\nproduct_daily = pd.concat(\n    product_daily_parts,\n    ignore_index=True\n)\n\n# Combine any duplicate date/product records\nproduct_daily = (\n    product_daily\n    .groupby(\n        [\"date\", \"item_nbr\"]\n    )[\"unit_sales\"]\n    .sum()\n    .reset_index()\n)\n\n\nprint(\"Rows extracted:\", len(product_daily))\n\ndisplay(\n    product_daily.head()\n)\n\n\n# ============================================================\n# 3. CREATE COMPLETE DAILY PRODUCT SERIES\n# ============================================================\n\nfull_dates = pd.date_range(\n    start=forecast_df[\"date\"].min(),\n    end=forecast_df[\"date\"].max(),\n    freq=\"D\"\n)\n\nproduct_series = (\n    product_daily\n    .pivot(\n        index=\"date\",\n        columns=\"item_nbr\",\n        values=\"unit_sales\"\n    )\n    .reindex(full_dates)\n)\n\n\nproduct_series.index.name = \"date\"\n\n\n# IMPORTANT:\n# Favorita does not explicitly store every zero-sales\n# product/date combination.\n#\n# Here we use zero for an unrecorded top-product/date\n# observation so that the product series can be modeled\n# at daily frequency.\n#\n# This is an assumption and should be documented in\n# the final report.\n\nproduct_series = product_series.fillna(0)\n\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"PRODUCT TIME SERIES\")\nprint(\"=\" * 70)\n\nprint(\"Shape:\", product_series.shape)\n\ndisplay(\n    product_series.head()\n)\n\n\n# ============================================================\n# 4. PRODUCT-LEVEL TEST PERIOD\n# ============================================================\n\nproduct_accuracy_results = []\n\nproduct_predictions = {}\n\n\nTEST_START = pd.Timestamp(\"2016-08-16\")\nTEST_END = pd.Timestamp(\"2017-08-15\")\n\n\n# ============================================================\n# 5. FORECAST EACH TOP PRODUCT\n# ============================================================\n\nfor item_id in top_product_ids:\n\n    print(\"\\n\" + \"=\" * 70)\n    print(\"PRODUCT:\", item_id)\n    print(\"=\" * 70)\n\n    product_ts = (\n        product_series[item_id]\n        .asfreq(\"D\")\n    )\n\n    product_train = product_ts.loc[\n        product_ts.index < TEST_START\n    ]\n\n    product_test = product_ts.loc[\n        (product_ts.index >= TEST_START) &\n        (product_ts.index <= TEST_END)\n    ]\n\n\n    print(\n        \"Training observations:\",\n        len(product_train)\n    )\n\n    print(\n        \"Test observations:\",\n        len(product_test)\n    )\n\n\n    # --------------------------------------------------------\n    # Seasonal Naive\n    # --------------------------------------------------------\n\n    combined_product = pd.concat([\n        product_train,\n        product_test\n    ])\n\n    seasonal_naive_pred = (\n        combined_product\n        .shift(7)\n        .loc[product_test.index]\n        .clip(lower=0)\n    )\n\n\n    # --------------------------------------------------------\n    # Tuned SARIMA\n    # --------------------------------------------------------\n\n    try:\n\n        product_model = SARIMAX(\n            product_train,\n            order=(2, 1, 2),\n            seasonal_order=(1, 0, 1, 7),\n            enforce_stationarity=False,\n            enforce_invertibility=False\n        )\n\n        product_result = product_model.fit(\n            disp=False\n        )\n\n        sarima_pred = (\n            product_result\n            .forecast(\n                steps=len(product_test)\n            )\n        )\n\n        sarima_pred.index = (\n            product_test.index\n        )\n\n        # Demand cannot be negative\n        sarima_pred = sarima_pred.clip(\n            lower=0\n        )\n\n\n        # ----------------------------------------------------\n        # Calculate metrics\n        # ----------------------------------------------------\n\n        sarima_metrics = calculate_metrics(\n            product_test,\n            sarima_pred,\n            f\"SARIMA - {item_id}\"\n        )\n\n        seasonal_metrics = calculate_metrics(\n            product_test,\n            seasonal_naive_pred,\n            f\"Seasonal Naive - {item_id}\"\n        )\n\n\n        # Add product metadata\n        product_family = (\n            top_5_products\n            .loc[\n                top_5_products[\"item_nbr\"] == item_id,\n                \"family\"\n            ]\n            .iloc[0]\n        )\n\n\n        for result in [\n            sarima_metrics,\n            seasonal_metrics\n        ]:\n\n            result[\"item_nbr\"] = item_id\n            result[\"family\"] = product_family\n\n\n        product_accuracy_results.extend([\n            sarima_metrics,\n            seasonal_metrics\n        ])\n\n\n        product_predictions[item_id] = {\n            \"actual\": product_test,\n            \"sarima\": sarima_pred,\n            \"seasonal_naive\":\n                seasonal_naive_pred\n        }\n\n\n        print(\"SARIMA MAE:\",\n              round(sarima_metrics[\"MAE\"], 2))\n\n        print(\"Seasonal Naive MAE:\",\n              round(seasonal_metrics[\"MAE\"], 2))\n\n\n    except Exception as e:\n\n        print(\n            \"SARIMA failed for product:\",\n            item_id\n        )\n\n        print(\n            \"Error:\",\n            str(e)[:200]\n        )\n\n\n# ============================================================\n# 6. PRODUCT ACCURACY TABLE\n# ============================================================\n\nproduct_accuracy = pd.DataFrame(\n    product_accuracy_results\n)\n\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"PRODUCT-LEVEL MODEL ACCURACY\")\nprint(\"=\" * 70)\n\ndisplay(\n    product_accuracy[\n        [\n            \"item_nbr\",\n            \"family\",\n            \"Model\",\n            \"MAE\",\n            \"RMSE\",\n            \"MAPE (%)\",\n            \"Directional Accuracy (%)\"\n        ]\n    ].round(3)\n)\n\n\n# ============================================================\n# 7. PRODUCT COMPARISON — SARIMA VS SEASONAL NAIVE\n# ============================================================\n\nproduct_comparison = (\n    product_accuracy\n    .pivot_table(\n        index=[\n            \"item_nbr\",\n            \"family\"\n        ],\n        columns=\"Model\",\n        values=[\n            \"MAE\",\n            \"RMSE\",\n            \"MAPE (%)\",\n            \"Directional Accuracy (%)\"\n        ]\n    )\n)\n\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"PRODUCT-LEVEL COMPARISON\")\nprint(\"=\" * 70)\n\ndisplay(\n    product_comparison.round(3)\n)\n\n\n# ============================================================\n# 8. ACTUAL VS PREDICTED — EACH PRODUCT\n# ============================================================\n\nfor item_id in top_product_ids:\n\n    if item_id not in product_predictions:\n        continue\n\n    result = product_predictions[item_id]\n\n    plt.figure(figsize=(16, 5))\n\n    plt.plot(\n        result[\"actual\"].index,\n        result[\"actual\"].values,\n        label=\"Actual\"\n    )\n\n    plt.plot(\n        result[\"sarima\"].index,\n        result[\"sarima\"].values,\n        label=\"SARIMA\"\n    )\n\n    plt.plot(\n        result[\"seasonal_naive\"].index,\n        result[\"seasonal_naive\"].values,\n        label=\"Seasonal Naive\"\n    )\n\n    family = (\n        top_5_products\n        .loc[\n            top_5_products[\"item_nbr\"] == item_id,\n            \"family\"\n        ]\n        .iloc[0]\n    )\n\n    plt.title(\n        f\"Product {item_id} ({family}) — \"\n        f\"Actual vs Forecast\"\n    )\n\n    plt.xlabel(\"Date\")\n    plt.ylabel(\"Unit Sales\")\n\n    plt.legend()\n\n    plt.xticks(rotation=45)\n\n    plt.tight_layout()\n    plt.show()\n\n\n# ============================================================\n# 9. FIND PRODUCTS WHERE SARIMA BEATS SEASONAL NAIVE\n# ============================================================\n\nsarima_product_mae = (\n    product_accuracy[\n        product_accuracy[\"Model\"]\n        .str.startswith(\"SARIMA\")\n    ][\n        [\n            \"item_nbr\",\n            \"MAE\"\n        ]\n    ]\n    .rename(\n        columns={\n            \"MAE\":\n                \"SARIMA_MAE\"\n        }\n    )\n)\n\nnaive_product_mae = (\n    product_accuracy[\n        product_accuracy[\"Model\"]\n        .str.startswith(\"Seasonal Naive\")\n    ][\n        [\n            \"item_nbr\",\n            \"MAE\"\n        ]\n    ]\n    .rename(\n        columns={\n            \"MAE\":\n                \"Seasonal_Naive_MAE\"\n        }\n    )\n)\n\nproduct_model_comparison = (\n    sarima_product_mae\n    .merge(\n        naive_product_mae,\n        on=\"item_nbr\",\n        how=\"inner\"\n    )\n)\n\nproduct_model_comparison[\n    \"SARIMA_better\"\n] = (\n    product_model_comparison[\n        \"SARIMA_MAE\"\n    ]\n    <\n    product_model_comparison[\n        \"Seasonal_Naive_MAE\"\n    ]\n)\n\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"PRODUCT MODEL WIN/LOSS\")\nprint(\"=\" * 70)\n\ndisplay(\n    product_model_comparison.round(3)\n)\n\n\nprint(\n    \"\\nSARIMA better for\",\n    product_model_comparison[\"SARIMA_better\"].sum(),\n    \"of\",\n    len(product_model_comparison),\n    \"products.\"\n)\n\n\n# ============================================================\n# 10. SAVE PRODUCT-LEVEL OUTPUTS\n# ============================================================\n\ntop_5_products.to_csv(\n    \"/kaggle/working/top_5_products.csv\",\n    index=False\n)\n\nproduct_daily.to_csv(\n    \"/kaggle/working/top_5_product_daily_sales.csv\",\n    index=False\n)\n\nproduct_accuracy.to_csv(\n    \"/kaggle/working/product_level_accuracy.csv\",\n    index=False\n)\n\nproduct_model_comparison.to_csv(\n    \"/kaggle/working/product_model_comparison.csv\",\n    index=False\n)\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"FILES SAVED\")\nprint(\"=\" * 70)\n\nprint(\"/kaggle/working/top_5_products.csv\")\nprint(\"/kaggle/working/top_5_product_daily_sales.csv\")\nprint(\"/kaggle/working/product_level_accuracy.csv\")\nprint(\"/kaggle/working/product_model_comparison.csv\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:33:09.081659Z","iopub.execute_input":"2026-09-17T10:33:09.083085Z","iopub.status.idle":"2026-09-17T10:35:40.784688Z","shell.execute_reply.started":"2026-09-17T10:33:09.083038Z","shell.execute_reply":"2026-09-17T10:35:40.783379Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# STEP 31 — FINAL PROJECT REPORT TABLES\n# ============================================================\n\nimport pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\n\n\n# ============================================================\n# 1. FINAL MODEL COMPARISON\n# ============================================================\n\nfinal_model_comparison = pd.concat(\n    [\n        baseline_results,\n        arima_results,\n        sarima_results,\n        exp_results,\n        prophet_results,\n        lstm_results\n    ],\n    ignore_index=True\n)\n\nfinal_model_comparison = (\n    final_model_comparison[\n        [\n            \"Model\",\n            \"MAE\",\n            \"RMSE\",\n            \"MAPE (%)\",\n            \"Directional Accuracy (%)\"\n        ]\n    ]\n    .sort_values(\"MAE\")\n    .reset_index(drop=True)\n)\n\n\nprint(\"=\" * 70)\nprint(\"FINAL MODEL COMPARISON — VALIDATION\")\nprint(\"=\" * 70)\n\ndisplay(\n    final_model_comparison.round(3)\n)\n\n\n# ============================================================\n# 2. FINAL TEST RESULTS\n# ============================================================\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"FINAL TEST RESULTS\")\nprint(\"=\" * 70)\n\ndisplay(\n    final_test_results.round(3)\n)\n\n\n# ============================================================\n# 3. PRODUCT-LEVEL ACCURACY\n# ============================================================\n\nproduct_accuracy_summary = (\n    product_accuracy[\n        [\n            \"item_nbr\",\n            \"family\",\n            \"Model\",\n            \"MAE\",\n            \"RMSE\",\n            \"MAPE (%)\",\n            \"Directional Accuracy (%)\"\n        ]\n    ]\n    .sort_values(\n        [\"item_nbr\", \"MAE\"]\n    )\n    .reset_index(drop=True)\n)\n\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"PRODUCT-LEVEL ACCURACY\")\nprint(\"=\" * 70)\n\ndisplay(\n    product_accuracy_summary.round(3)\n)\n\n\n# ============================================================\n# 4. PRODUCT MODEL WIN COUNT\n# ============================================================\n\nproduct_win_summary = pd.DataFrame({\n    \"Comparison\": [\n        \"SARIMA better than Seasonal Naive\"\n    ],\n    \"Products Won\": [\n        product_model_comparison[\n            \"SARIMA_better\"\n        ].sum()\n    ],\n    \"Products Evaluated\": [\n        len(product_model_comparison)\n    ]\n})\n\nproduct_win_summary[\"Win Rate (%)\"] = (\n    product_win_summary[\"Products Won\"]\n    /\n    product_win_summary[\"Products Evaluated\"]\n    * 100\n)\n\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"PRODUCT-LEVEL MODEL WIN RATE\")\nprint(\"=\" * 70)\n\ndisplay(\n    product_win_summary.round(2)\n)\n\n\n# ============================================================\n# 5. FORECAST SUMMARY\n# ============================================================\n\nforecast_summary = pd.DataFrame({\n\n    \"Horizon\": [\n        \"6 Months\",\n        \"12 Months\"\n    ],\n\n    \"Start Date\": [\n        forecast_6m.index.min(),\n        forecast_12m.index.min()\n    ],\n\n    \"End Date\": [\n        forecast_6m.index.max(),\n        forecast_12m.index.max()\n    ],\n\n    \"Total Forecast Demand\": [\n        forecast_6m[\"forecast\"].sum(),\n        forecast_12m[\"forecast\"].sum()\n    ],\n\n    \"Average Daily Demand\": [\n        forecast_6m[\"forecast\"].mean(),\n        forecast_12m[\"forecast\"].mean()\n    ]\n})\n\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"FINAL FORECAST SUMMARY\")\nprint(\"=\" * 70)\n\ndisplay(\n    forecast_summary.round(2)\n)\n\n\n# ============================================================\n# 6. INVENTORY SUMMARY\n# ============================================================\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"INVENTORY PLANNING SUMMARY\")\nprint(\"=\" * 70)\n\ndisplay(\n    inventory_planning.round(2)\n)\n\n\n# ============================================================\n# 7. TOP PRODUCTS\n# ============================================================\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"TOP 5 PRODUCTS BY HISTORICAL DEMAND\")\nprint(\"=\" * 70)\n\ndisplay(\n    top_5_products.round(2)\n)\n\n\n# ============================================================\n# 8. BEST VALIDATION MODEL\n# ============================================================\n\nbest_validation_model = (\n    final_model_comparison\n    .iloc[0][\"Model\"]\n)\n\nbest_validation_mae = (\n    final_model_comparison\n    .iloc[0][\"MAE\"]\n)\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"BEST VALIDATION MODEL BY MAE\")\nprint(\"=\" * 70)\n\nprint(\n    \"Model:\",\n    best_validation_model\n)\n\nprint(\n    \"MAE:\",\n    round(best_validation_mae, 3)\n)\n\n\n# ============================================================\n# 9. FINAL TEST MODEL COMPARISON\n# ============================================================\n\nbest_test_by_mae = (\n    final_test_results\n    .sort_values(\"MAE\")\n    .iloc[0]\n)\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"BEST TEST MODEL BY MAE\")\nprint(\"=\" * 70)\n\nprint(\n    \"Model:\",\n    best_test_by_mae[\"Model\"]\n)\n\nprint(\n    \"MAE:\",\n    round(best_test_by_mae[\"MAE\"], 3)\n)\n\nprint(\n    \"RMSE:\",\n    round(best_test_by_mae[\"RMSE\"], 3)\n)\n\nprint(\n    \"MAPE:\",\n    round(best_test_by_mae[\"MAPE (%)\"], 3)\n)\n\nprint(\n    \"Directional Accuracy:\",\n    round(\n        best_test_by_mae[\n            \"Directional Accuracy (%)\"\n        ],\n        3\n    )\n)\n\n\n# ============================================================\n# 10. FINAL FORECAST PLOT\n# ============================================================\n\nplt.figure(figsize=(18, 7))\n\nrecent_history = final_series.tail(180)\n\nplt.plot(\n    recent_history.index,\n    recent_history.values,\n    label=\"Historical\"\n)\n\nplt.plot(\n    forecast_output.index,\n    forecast_output[\"forecast\"],\n    label=\"12-Month Forecast\"\n)\n\nplt.fill_between(\n    forecast_output.index,\n    forecast_output[\"lower_95_inventory\"],\n    forecast_output[\"upper_95_inventory\"],\n    alpha=0.2,\n    label=\"95% Confidence Interval\"\n)\n\nplt.axvline(\n    final_series.index.max(),\n    linestyle=\"--\",\n    label=\"Forecast Start\"\n)\n\nplt.title(\n    \"Final Demand Forecast — Historical + 12 Months\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\nplt.xticks(rotation=45)\n\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 11. MODEL MAE COMPARISON\n# ============================================================\n\nplt.figure(figsize=(12, 6))\n\nplt.bar(\n    final_model_comparison[\"Model\"],\n    final_model_comparison[\"MAE\"]\n)\n\nplt.title(\n    \"Forecasting Model Comparison — MAE\"\n)\n\nplt.xlabel(\"Model\")\nplt.ylabel(\"MAE\")\n\nplt.xticks(\n    rotation=45,\n    ha=\"right\"\n)\n\nplt.tight_layout()\nplt.show()\n\n\n# ============================================================\n# 12. SAVE FINAL REPORT TABLES\n# ============================================================\n\nfinal_model_comparison.to_csv(\n    \"/kaggle/working/final_model_comparison.csv\",\n    index=False\n)\n\nfinal_test_results.to_csv(\n    \"/kaggle/working/final_test_model_results.csv\",\n    index=False\n)\n\nproduct_accuracy_summary.to_csv(\n    \"/kaggle/working/final_product_accuracy.csv\",\n    index=False\n)\n\nproduct_win_summary.to_csv(\n    \"/kaggle/working/product_model_win_summary.csv\",\n    index=False\n)\n\nforecast_summary.to_csv(\n    \"/kaggle/working/final_forecast_summary.csv\",\n    index=False\n)\n\ninventory_planning.to_csv(\n    \"/kaggle/working/final_inventory_planning.csv\",\n    index=False\n)\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"FINAL REPORT FILES SAVED\")\nprint(\"=\" * 70)\n\nprint(\"/kaggle/working/final_model_comparison.csv\")\nprint(\"/kaggle/working/final_test_model_results.csv\")\nprint(\"/kaggle/working/final_product_accuracy.csv\")\nprint(\"/kaggle/working/product_model_win_summary.csv\")\nprint(\"/kaggle/working/final_forecast_summary.csv\")\nprint(\"/kaggle/working/final_inventory_planning.csv\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:36:57.521842Z","iopub.execute_input":"2026-09-17T10:36:57.522325Z","iopub.status.idle":"2026-09-17T10:36:58.321614Z","shell.execute_reply.started":"2026-09-17T10:36:57.522289Z","shell.execute_reply":"2026-09-17T10:36:58.320262Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# STEP 33 — FINAL DELIVERABLES REPORT\n# ============================================================\n#\n# Creates ONE Markdown report containing:\n#\n# 1. Project overview\n# 2. Dataset and methodology\n# 3. EDA findings\n# 4. Stationarity + ACF/PACF\n# 5. Model comparison\n# 6. Hyperparameter tuning\n# 7. Backtesting\n# 8. Final test evaluation\n# 9. Actual vs predicted plot placeholders\n# 10. Product-level accuracy metrics\n# 11. Confidence interval plot placeholders\n# 12. 6/12-month forecasts\n# 13. Inventory recommendations\n# 14. Business recommendations\n# 15. Limitations\n# 16. Conclusion\n#\n# Output:\n# /kaggle/working/Favorita_Demand_Forecasting_Report.md\n# ============================================================\n\nfrom pathlib import Path\nimport pandas as pd\nimport numpy as np\n\n\n# ============================================================\n# 1. HELPER FUNCTION\n# ============================================================\n\ndef md_table(df, index=False, round_digits=2):\n    \"\"\"\n    Convert DataFrame to clean Markdown table.\n    \"\"\"\n    temp = df.copy()\n\n    # Round numeric columns\n    for col in temp.select_dtypes(\n        include=[\"float64\", \"float32\", \"int64\", \"int32\"]\n    ).columns:\n\n        temp[col] = temp[col].round(round_digits)\n\n    try:\n        return temp.to_markdown(\n            index=index\n        )\n    except Exception:\n        return temp.to_string(\n            index=index\n        )\n\n\n# ============================================================\n# 2. REPORT PATH\n# ============================================================\n\nreport_path = Path(\n    \"/kaggle/working/\"\n    \"Favorita_Demand_Forecasting_Report.md\"\n)\n\n\n# ============================================================\n# 3. SAFELY PREPARE AVAILABLE RESULTS\n# ============================================================\n\n# Validation model comparison\nvalidation_results = final_model_comparison.copy()\n\n# Final test results\ntest_results = final_test_results.copy()\n\n# Product-level accuracy\nproduct_results = product_accuracy_summary.copy()\n\n# Forecast summary\nforecast_results = forecast_summary.copy()\n\n# Inventory planning\ninventory_results = inventory_planning.copy()\n\n# Top products\ntop_products = top_5_products.copy()\n\n# Monthly forecasts\nmonthly_results = monthly_forecast.copy()\n\n\n# ============================================================\n# 4. EXTRACT IMPORTANT VALUES\n# ============================================================\n\ntry:\n\n    tuned_sarima_validation = best_sarima_metrics_df.iloc[0]\n\nexcept Exception:\n\n    tuned_sarima_validation = sarima_results.iloc[0]\n\n\ntry:\n\n    test_sarima = test_results[\n        test_results[\"Model\"] == \"Tuned SARIMA\"\n    ].iloc[0]\n\nexcept Exception:\n\n    test_sarima = test_results.iloc[0]\n\n\ntry:\n\n    test_seasonal_naive = test_results[\n        test_results[\"Model\"] == \"Seasonal Naive\"\n    ].iloc[0]\n\nexcept Exception:\n\n    test_seasonal_naive = None\n\n\n# Forecast quantities\nsix_month_total = forecast_6m[\"forecast\"].sum()\nsix_month_average = forecast_6m[\"forecast\"].mean()\n\ntwelve_month_total = forecast_12m[\"forecast\"].sum()\ntwelve_month_average = forecast_12m[\"forecast\"].mean()\n\n\n# Highest forecast month\nhighest_forecast_month = monthly_results.loc[\n    monthly_results[\"forecast_units\"].idxmax()\n]\n\n\n# Lowest forecast month\nlowest_forecast_month = monthly_results.loc[\n    monthly_results[\"forecast_units\"].idxmin()\n]\n\n\n# Product-level SARIMA win count\ntry:\n\n    product_sarima_wins = int(\n        product_model_comparison[\n            \"SARIMA_better\"\n        ].sum()\n    )\n\n    product_total_evaluated = int(\n        len(product_model_comparison)\n    )\n\n    product_win_rate = (\n        product_sarima_wins\n        /\n        product_total_evaluated\n        *\n        100\n    )\n\nexcept Exception:\n\n    product_sarima_wins = None\n    product_total_evaluated = None\n    product_win_rate = None\n\n\n# ============================================================\n# 5. BUILD REPORT\n# ============================================================\n\nreport = f\"\"\"\n# Favorita Grocery Sales Demand Forecasting\n\n## Model Comparison, Forecasting Results, Accuracy Analysis, and Inventory Planning\n\n---\n\n# 1. Executive Summary\n\n### Objective\n\nThe objective of this project is to forecast future product demand\nfor approximately 6–12 months using the Favorita Grocery Sales\nForecasting dataset.\n\nThe forecasts are intended to support:\n\n- inventory planning\n- supply-chain planning\n- replenishment decisions\n- seasonal procurement\n- demand-risk assessment\n\nThe main aggregate forecasting series was constructed by\naggregating sales by date across stores and products.\n\n### Historical Period\n\n**2013-01-01 to 2017-08-15**\n\n### Main Forecast Horizon\n\n**6 months:** {forecast_6m.index.min().date()} to {forecast_6m.index.max().date()}\n\n**12 months:** {forecast_12m.index.min().date()} to {forecast_12m.index.max().date()}\n\n### Final Forecast Model\n\nTuned seasonal ARIMA:\n\n**SARIMA(2,1,2)(1,0,1,7)**\n\n### Important Evaluation Finding\n\nThe tuned SARIMA model produced the lowest MAE on the final\ntest period, while the 7-day Seasonal Naive model produced lower\nMAPE and higher directional accuracy.\n\nTherefore, no single model dominates every evaluation metric.\n\n---\n\n# 2. Dataset Overview\n\n## Dataset\n\nFavorita Grocery Sales Forecasting dataset from Kaggle.\n\n### Main datasets used\n\n- `train.csv`\n- `items.csv`\n- `stores.csv`\n- `holidays_events.csv`\n\nAdditional Favorita datasets were available for extended feature\nengineering.\n\n## Dataset scale\n\n### Product catalog\n\n- 4,100 products\n- 33 product families\n- 337 product classes\n\n### Stores\n\n- 54 stores\n- 22 cities\n- 16 states\n- 5 store types\n- 17 store clusters\n\n### Historical sales\n\nThe extracted training CSV was approximately 5 GB, so chunk-based\nprocessing was used instead of loading the complete file into memory.\n\n---\n\n# 3. Data Preparation\n\nThe original sales data was recorded approximately at:\n\n**Date × Store × Product**\n\nlevel.\n\nThe aggregate forecasting dataset was created as:\n\n**Date → Total Unit Sales**\n\nThis produced a daily time series covering the historical period.\n\n## Missing calendar dates\n\nFour dates were absent from the raw aggregate sales series:\n\n- 2013-12-25\n- 2014-12-25\n- 2015-12-25\n- 2016-12-25\n\nThese corresponded to Christmas (`Navidad`) national holidays in\nthe holiday dataset.\n\nFor the aggregate forecasting series, these dates were represented\nas zero recorded demand and retained with a `structural_missing`\nindicator.\n\n## Negative sales\n\nThe raw data contained negative `unit_sales` observations,\nconsistent with returns in the competition data.\n\nNumber of negative raw observations identified:\n\n**7,795**\n\nThese were not blindly removed before aggregation.\n\n---\n\n# 4. Exploratory Data Analysis\n\n## Major EDA Findings\n\n### Long-term trend\n\nAverage daily demand increased substantially across the historical\nperiod.\n\n### Monthly seasonality\n\nDecember showed the highest average daily sales among the months.\n\n### Weekly seasonality\n\nDemand varied considerably by day of week.\n\nThe highest average daily demand occurred on:\n\n**Sunday**\n\nThe lowest occurred on:\n\n**Thursday**\n\n### Holiday effect\n\nHoliday and non-holiday demand showed different average sales\nlevels.\n\n### Outliers\n\nFive aggregate daily observations were flagged using the IQR\nmethod.\n\nThese observations were retained because unusually high demand\ncan represent genuine business behavior rather than bad data.\n\n---\n\n# 5. EDA Visualizations\n\n## Daily Sales Trend\n\n[image - Daily Total Sales Over Time]\n\n## Daily Sales Distribution\n\n[image - Distribution of Daily Total Sales]\n\n## Weekly Sales Trend\n\n[image - Weekly Total Sales Over Time]\n\n## Monthly Sales Trend\n\n[image - Monthly Total Sales Over Time]\n\n## Average Sales by Month\n\n[image - Average Daily Sales by Month]\n\n## Average Sales by Day of Week\n\n[image - Average Daily Sales by Day of Week]\n\n## Year × Month Seasonality\n\n[image - Year by Month Sales Heatmap]\n\n## Holiday vs Non-Holiday Demand\n\n[image - Holiday vs Non-Holiday Sales]\n\n---\n\n# 6. Time-Based Feature Engineering\n\nThe final feature dataset included:\n\n- year\n- month\n- quarter\n- day\n- day of week\n- day name\n- weekend indicator\n- holiday indicator\n- holiday count\n- holiday type\n- holiday description\n- holiday locale\n\nThe dataset was saved as:\n\n`daily_sales_features.csv`\n\n---\n\n# 7. Stationarity Analysis\n\nThe Augmented Dickey-Fuller test was performed.\n\n## Original Series\n\nADF Statistic:\n\n**-2.5721**\n\np-value:\n\n**0.09891**\n\nConclusion:\n\nThe original series was treated as non-stationary at the 5%\nsignificance level.\n\n## First Difference\n\nADF Statistic:\n\n**-11.3137**\n\np-value:\n\n**1.22e-20**\n\nConclusion:\n\nThe first-differenced series was stationary.\n\n## Seven-Day Difference\n\nADF Statistic:\n\n**-12.7575**\n\np-value:\n\n**8.28e-24**\n\nConclusion:\n\nThe 7-day differenced series was stationary.\n\n---\n\n# 8. ACF and PACF Analysis\n\nThe first-differenced ACF showed strong autocorrelation at lag 7.\n\n### Important result\n\n**ACF at lag 7 ≈ 0.646**\n\nThis provided evidence of weekly seasonality.\n\nThe forecasting models therefore used:\n\n**Seasonal period = 7 days**\n\n## ACF Plot\n\n[image - ACF Plot]\n\n## PACF Plot\n\n[image - PACF Plot]\n\n## Seasonal ACF Plot\n\n[image - Seasonal ACF Plot]\n\n## Seasonal PACF Plot\n\n[image - Seasonal PACF Plot]\n\n---\n\n# 9. Train / Validation / Test Strategy\n\nThe data was split chronologically.\n\n## Training\n\n**2013-01-01 to 2016-02-15**\n\n1,141 observations.\n\n## Validation\n\n**2016-02-16 to 2016-08-15**\n\n182 observations.\n\n## Test\n\n**2016-08-16 to 2017-08-15**\n\n365 observations.\n\nNo random shuffling was used.\n\nThis preserves the time-dependent nature of forecasting and avoids\nfuture-information leakage.\n\n---\n\n# 10. Baseline Models\n\nTwo simple baselines were established:\n\n### Naive\n\nThe next value is based on the latest observed value.\n\n### Seasonal Naive\n\nThe forecast uses the value from 7 days earlier.\n\nThe Seasonal Naive baseline was particularly important because\nweekly seasonality was detected in the EDA and ACF analysis.\n\n---\n\n# 11. Forecasting Models\n\nThe following models were implemented:\n\n1. Naive\n2. Seasonal Naive\n3. ARIMA\n4. SARIMA\n5. Exponential Smoothing\n6. Prophet\n7. LSTM\n\nEvaluation metrics:\n\n- MAE\n- RMSE\n- MAPE\n- Directional Accuracy\n\n---\n\n# 12. Model Comparison — Validation\n\n{md_table(validation_results, index=False, round_digits=3)}\n\n### Validation Interpretation\n\nThe initial SARIMA model:\n\n**SARIMA(1,1,1)(1,0,1,7)**\n\nproduced the lowest validation MAE:\n\n**{sarima_results[\"MAE\"].iloc[0]:,.3f}**\n\nThe Exponential Smoothing model was also competitive.\n\nThe LSTM model produced higher error on this aggregate daily series\nthan the statistical seasonal models.\n\n---\n\n# 13. Hyperparameter Tuning\n\nA SARIMA grid search containing:\n\n**36 parameter combinations**\n\nwas performed.\n\nThe selected tuned model was:\n\n**SARIMA(2,1,2)(1,0,1,7)**\n\n## Tuned Validation Metrics\n\n| Metric | Result |\n|---|---:|\n| MAE | {tuned_sarima_validation[\"MAE\"]:,.3f} |\n| RMSE | {tuned_sarima_validation[\"RMSE\"]:,.3f} |\n| MAPE | {tuned_sarima_validation[\"MAPE (%)\"]:.3f}% |\n| Directional Accuracy | {tuned_sarima_validation[\"Directional Accuracy (%)\"]:.3f}% |\n\nThe tuning improvement in MAE was approximately:\n\n**0.3%**\n\n---\n\n# 14. Backtesting\n\nExpanding-window backtesting was performed using:\n\n- 90-day forecast horizon\n- 4 historical folds\n- tuned SARIMA\n- 7-day Seasonal Naive\n\n## Pooled Backtesting Results\n\n{md_table(\n    pooled_results[\n        [\n            \"Model\",\n            \"MAE\",\n            \"RMSE\",\n            \"MAPE (%)\",\n            \"Directional Accuracy (%)\"\n        ]\n    ],\n    index=False,\n    round_digits=3\n)}\n\n### Backtesting Interpretation\n\nThe Seasonal Naive model was more consistent across the four\nhistorical folds.\n\nTuned SARIMA performed better on:\n\n**2 of 4 folds**\n\nThis demonstrates that model performance changes depending on the\nhistorical period.\n\n## Backtesting MAE\n\n[image - Backtesting MAE by Fold]\n\n## Backtesting Example\n\n[image - Backtest Fold 1 Actual vs Forecast]\n\n---\n\n# 15. Final Test Evaluation\n\nThe final test period was:\n\n**2016-08-16 to 2017-08-15**\n\n{md_table(test_results, index=False, round_digits=3)}\n\n## Test Interpretation\n\n### Tuned SARIMA\n\n- MAE: {test_sarima[\"MAE\"]:,.3f}\n- RMSE: {test_sarima[\"RMSE\"]:,.3f}\n- MAPE: {test_sarima[\"MAPE (%)\"]:.3f}%\n- Directional Accuracy: {test_sarima[\"Directional Accuracy (%)\"]:.3f}%\n\n### Seasonal Naive\n\n\"\"\"\n\nif test_seasonal_naive is not None:\n\n    report += f\"\"\"\n- MAE: {test_seasonal_naive[\"MAE\"]:,.3f}\n- RMSE: {test_seasonal_naive[\"RMSE\"]:,.3f}\n- MAPE: {test_seasonal_naive[\"MAPE (%)\"]:.3f}%\n- Directional Accuracy: {test_seasonal_naive[\"Directional Accuracy (%)\"]:.3f}%\n\"\"\"\n\nreport += f\"\"\"\nThe tuned SARIMA model achieved a lower MAE on the final test\nperiod, while Seasonal Naive achieved lower MAPE and higher\ndirectional accuracy.\n\nThis difference should be explicitly considered when interpreting\nforecast accuracy for business use.\n\n---\n\n# 16. Final Actual vs Predicted Visualizations\n\n## Validation — Seasonal Naive\n\n[image - Actual vs Seasonal Naive — Validation]\n\n## Validation — ARIMA\n\n[image - Actual vs ARIMA — Validation]\n\n## Validation — SARIMA\n\n[image - Actual vs SARIMA — Validation]\n\n## Validation — Exponential Smoothing\n\n[image - Actual vs Exponential Smoothing — Validation]\n\n## Validation — Prophet\n\n[image - Actual vs Prophet — Validation]\n\n## Validation — LSTM\n\n[image - Actual vs LSTM — Validation]\n\n## Validation — Main Model Comparison\n\n[image - Actual vs SARIMA vs Exponential Smoothing vs Prophet vs LSTM]\n\n## Final Test — Actual vs Tuned SARIMA\n\n[image - Final Test Actual vs Tuned SARIMA]\n\n## Final Test — Actual vs Seasonal Naive\n\n[image - Final Test Actual vs Seasonal Naive]\n\n---\n\n# 17. Product-Level Accuracy\n\nThe product-level analysis evaluated five high-demand products.\n\n{md_table(\n    product_results[\n        [\n            \"item_nbr\",\n            \"family\",\n            \"Model\",\n            \"MAE\",\n            \"RMSE\",\n            \"MAPE (%)\",\n            \"Directional Accuracy (%)\"\n        ]\n    ],\n    index=False,\n    round_digits=3\n)}\n\n## Product-Level Model Comparison\n\n[image - Product-Level Model Accuracy Comparison]\n\n### Product-Level Findings\n\nSARIMA achieved lower MAE than Seasonal Naive for:\n\n**{product_sarima_wins} of {product_total_evaluated} products**\n\nWin rate:\n\n**{product_win_rate:.1f}%**\n\nThe results show that the best forecasting method can vary by\nproduct.\n\n---\n\n# 18. Product-Level Actual vs Predicted Plots\n\n## Product 1\n\n[image - Product 1 Actual vs Predicted]\n\n## Product 2\n\n[image - Product 2 Actual vs Predicted]\n\n## Product 3\n\n[image - Product 3 Actual vs Predicted]\n\n## Product 4\n\n[image - Product 4 Actual vs Predicted]\n\n## Product 5\n\n[image - Product 5 Actual vs Predicted]\n\n---\n\n# 19. Top Products by Historical Demand\n\n{md_table(\n    top_products,\n    index=False,\n    round_digits=2\n)}\n\n---\n\n# 20. Final 6-Month Forecast\n\nForecast period:\n\n**{forecast_6m.index.min().date()} to {forecast_6m.index.max().date()}**\n\nTotal forecast demand:\n\n**{six_month_total:,.2f} units**\n\nAverage daily forecast:\n\n**{six_month_average:,.2f} units**\n\n## 6-Month Forecast Visualization\n\n[image - 6-Month Demand Forecast with 95% Confidence Interval]\n\n---\n\n# 21. Final 12-Month Forecast\n\nForecast period:\n\n**{forecast_12m.index.min().date()} to {forecast_12m.index.max().date()}**\n\nTotal forecast demand:\n\n**{twelve_month_total:,.2f} units**\n\nAverage daily forecast:\n\n**{twelve_month_average:,.2f} units**\n\n## 12-Month Forecast Visualization\n\n[image - 12-Month Demand Forecast with 95% Confidence Interval]\n\n---\n\n# 22. Monthly Forecast\n\n{md_table(\n    monthly_results[\n        [\n            \"date\",\n            \"forecast_units\",\n            \"average_daily_forecast\",\n            \"lower_95_units\",\n            \"upper_95_units\"\n        ]\n    ],\n    index=False,\n    round_digits=2\n)}\n\n## Monthly Forecast Plot\n\n[image - Forecast Monthly Demand]\n\n---\n\n# 23. Forecast Confidence Intervals\n\nThe final SARIMA forecast includes 95% confidence intervals.\n\nThe intervals become wider as the forecast horizon increases,\nindicating increased uncertainty for longer-range forecasts.\n\n## Forecast with Confidence Interval\n\n[image - Historical Sales and 12-Month Forecast]\n\n## Six-Month Confidence Interval\n\n[image - 6-Month Forecast Confidence Interval]\n\n## Forecast Uncertainty\n\n[image - Forecast Uncertainty Over 12 Months]\n\n---\n\n# 24. Highest and Lowest Forecast-Demand Months\n\n### Highest Forecast Month\n\nDate:\n\n**{highest_forecast_month[\"date\"]}**\n\nForecast demand:\n\n**{highest_forecast_month[\"forecast_units\"]:,.2f} units**\n\n### Lowest Forecast Month\n\nDate:\n\n**{lowest_forecast_month[\"date\"]}**\n\nForecast demand:\n\n**{lowest_forecast_month[\"forecast_units\"]:,.2f} units**\n\n---\n\n# 25. Inventory Planning\n\n{md_table(\n    inventory_results,\n    index=False,\n    round_digits=2\n)}\n\n## Inventory Planning Visualization\n\n[image - Forecast Demand and Uncertainty for Inventory Planning]\n\n---\n\n# 26. Business Recommendations\n\n## 1. Demand Planning\n\nUse the forecasting model as the baseline demand plan and\nregularly update it as new sales observations become available.\n\n## 2. Weekly Replenishment\n\nThe strong weekly seasonal pattern supports replenishment\nplanning based on day-of-week demand differences.\n\n## 3. Seasonal Procurement\n\nProcurement and warehouse capacity should be reviewed ahead of\nhigher-demand forecast periods.\n\n## 4. Safety Stock\n\nThe upper forecast confidence bound can be used as a conservative\nplanning scenario when estimating inventory buffers.\n\nIt should not be treated as guaranteed future demand.\n\n## 5. Short-Term Replenishment\n\nThe 7-day and 14-day demand estimates can support short-term\ninventory planning when supplier lead times are relatively short.\n\n## 6. Longer Lead Times\n\nFor approximately 30-day supplier lead times, use the 30-day\nforecast together with the uncertainty buffer.\n\n## 7. Long-Range Planning\n\nThe uncertainty interval becomes wider as the forecast horizon\nincreases.\n\nTherefore, long-range decisions should be periodically\nre-forecast rather than locked in for the full 12 months.\n\n## 8. Product-Level Planning\n\nOperational inventory decisions should ultimately use\nproduct-store-level forecasts rather than only aggregate demand.\n\n## 9. Model Selection\n\nThe final analysis indicates that model selection should consider\nmultiple business metrics.\n\nMAE, MAPE, RMSE, and directional accuracy can tell different\nstories.\n\n---\n\n# 27. Project Limitations\n\n### Aggregate Forecasting\n\nThe primary SARIMA experiment forecasts total aggregate demand,\nnot every product-store combination independently.\n\n### Product-Level Coverage\n\nProduct-level evaluation was performed on five high-demand\nproducts rather than all 4,100 products.\n\n### External Variables\n\nThe final SARIMA model did not directly incorporate every available\nexternal variable such as promotion status, oil prices, store\ncharacteristics, and detailed holiday effects.\n\n### Product-Level Missingness Assumption\n\nWhen creating complete daily product-level series, unrecorded\nproduct/date observations were represented as zero sales.\n\nThis is a modeling assumption and should be documented when\ninterpreting those results.\n\n### Long-Term Forecast Uncertainty\n\nConfidence intervals widen substantially farther into the future.\n\nTherefore, the 12-month forecast should be interpreted as a\nplanning estimate rather than a precise prediction.\n\n---\n\n# 28. Recommended Future Improvements\n\nFuture versions of the project can improve accuracy by:\n\n1. Forecasting at Store × Product level.\n2. Adding promotion features.\n3. Adding oil-price variables.\n4. Incorporating detailed holiday/event variables.\n5. Using lag and rolling-demand features.\n6. Testing gradient boosting models.\n7. Using multivariate LSTM models.\n8. Testing LightGBM/XGBoost-style demand models.\n9. Performing hierarchical forecasting across product and store\n   levels.\n10. Automating periodic reforecasting.\n\n---\n\n# 29. Final Conclusion\n\nThis project developed a complete time-series forecasting pipeline\nfor the Favorita grocery sales dataset.\n\nThe workflow included:\n\n**Data ingestion**\n\n→ **Data cleaning**\n\n→ **Exploratory Data Analysis**\n\n→ **Trend and seasonality analysis**\n\n→ **Time-based feature engineering**\n\n→ **Stationarity testing**\n\n→ **Differencing**\n\n→ **ACF/PACF analysis**\n\n→ **Baseline forecasting**\n\n→ **ARIMA**\n\n→ **SARIMA**\n\n→ **Exponential Smoothing**\n\n→ **Prophet**\n\n→ **LSTM**\n\n→ **Hyperparameter tuning**\n\n→ **Historical backtesting**\n\n→ **Final test evaluation**\n\n→ **6-month forecasting**\n\n→ **12-month forecasting**\n\n→ **Confidence intervals**\n\n→ **Product-level accuracy**\n\n→ **Inventory planning**\n\nThe aggregate analysis identified strong weekly seasonality and\nsubstantial changes in demand over time.\n\nThe tuned SARIMA model achieved:\n\n- Test MAE: **{test_sarima[\"MAE\"]:,.2f}**\n- Test RMSE: **{test_sarima[\"RMSE\"]:,.2f}**\n- Test MAPE: **{test_sarima[\"MAPE (%)\"]:.2f}%**\n- Test Directional Accuracy: **{test_sarima[\"Directional Accuracy (%)\"]:.2f}%**\n\nThe final forecasting period covers approximately 6–12 months\nafter the end of the historical data.\n\nThese forecasts and their uncertainty intervals can support\naggregate demand planning, inventory buffering, replenishment\nplanning, and seasonal supply-chain preparation.\n\n---\n\n# 30. Deliverables Checklist\n\n| Deliverable | Status |\n|---|---|\n| Jupyter Notebook | Completed — export notebook separately |\n| Model Comparison Report | Included in this Markdown report |\n| Actual vs Predicted Visualizations | Placeholders included |\n| Accuracy Metrics by Product | Included |\n| Business Recommendations | Included |\n| Confidence Interval Plots | Placeholders included |\n| 6-Month Forecast | Included |\n| 12-Month Forecast | Included |\n| Backtesting | Included |\n| Hyperparameter Tuning | Included |\n| Stationarity Analysis | Included |\n| ACF/PACF Analysis | Included |\n\n---\n\n# 31. Generated Project Files\n\nThe notebook/exported project should include the following outputs:\n\n- `daily_sales_features.csv`\n- `final_forecasting_dataset.csv`\n- `train_timeseries.csv`\n- `validation_timeseries.csv`\n- `test_timeseries.csv`\n- `baseline_results.csv`\n- `arima_results.csv`\n- `sarima_results.csv`\n- `exponential_smoothing_results.csv`\n- `prophet_results.csv`\n- `lstm_results.csv`\n- `sarima_tuning_results.csv`\n- `backtest_results.csv`\n- `pooled_backtest_results.csv`\n- `final_test_results.csv`\n- `top_5_products.csv`\n- `product_level_accuracy.csv`\n- `final_6_month_forecast.csv`\n- `final_12_month_forecast.csv`\n- `monthly_forecast_summary.csv`\n- `final_inventory_planning.csv`\n\n---\n\n# END OF REPORT\n\"\"\"\n\n\n# ============================================================\n# 6. WRITE MARKDOWN FILE\n# ============================================================\n\nreport_path.write_text(\n    report,\n    encoding=\"utf-8\"\n)\n\n\n# ============================================================\n# 7. FINAL OUTPUT\n# ============================================================\n\nprint(\"=\" * 80)\nprint(\"FINAL DELIVERABLES REPORT CREATED\")\nprint(\"=\" * 80)\n\nprint(f\"\\nFile:\")\nprint(report_path)\n\nprint(\"\\nThe Markdown report contains:\")\nprint(\"✓ Project overview\")\nprint(\"✓ Data preparation\")\nprint(\"✓ EDA\")\nprint(\"✓ Stationarity\")\nprint(\"✓ ACF/PACF\")\nprint(\"✓ Model comparison\")\nprint(\"✓ Hyperparameter tuning\")\nprint(\"✓ Backtesting\")\nprint(\"✓ Final test evaluation\")\nprint(\"✓ Actual vs predicted image placeholders\")\nprint(\"✓ Product-level accuracy\")\nprint(\"✓ Product-level image placeholders\")\nprint(\"✓ 6-month forecast\")\nprint(\"✓ 12-month forecast\")\nprint(\"✓ Confidence interval image placeholders\")\nprint(\"✓ Inventory planning\")\nprint(\"✓ Business recommendations\")\nprint(\"✓ Limitations\")\nprint(\"✓ Conclusion\")\nprint(\"✓ Deliverables checklist\")\n\nprint(\"\\nUse the placeholders such as:\")\nprint(\"[image - 12-Month Demand Forecast with 95% Confidence Interval]\")\nprint(\"to paste the corresponding notebook plots below them.\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:44:08.461099Z","iopub.execute_input":"2026-09-17T10:44:08.465323Z","iopub.status.idle":"2026-09-17T10:44:08.637883Z","shell.execute_reply.started":"2026-09-17T10:44:08.465164Z","shell.execute_reply":"2026-09-17T10:44:08.636396Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# STEP 34 — GENERATE ALL REPORT IMAGES SEQUENTIALLY\n# ============================================================\n\nimport os\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n\nfrom statsmodels.graphics.tsaplots import plot_acf, plot_pacf\nfrom statsmodels.tsa.stattools import acf, pacf\n\n\n# ============================================================\n# 1. IMAGE OUTPUT DIRECTORY\n# ============================================================\n\nimage_dir = \"/kaggle/working/final_report_images\"\n\nos.makedirs(\n    image_dir,\n    exist_ok=True\n)\n\nprint(\"=\" * 80)\nprint(\"FINAL REPORT IMAGE GENERATION\")\nprint(\"=\" * 80)\nprint(f\"Images will be saved to:\\n{image_dir}\\n\")\n\n\n# ============================================================\n# 2. HELPER FUNCTION\n# ============================================================\n\nimage_counter = 0\n\n\ndef save_report_image(\n    title,\n    filename\n):\n    \"\"\"\n    Save current matplotlib figure with a sequential number.\n    \"\"\"\n\n    global image_counter\n\n    image_counter += 1\n\n    numbered_filename = (\n        f\"{image_counter:02d}_{filename}\"\n    )\n\n    output_path = os.path.join(\n        image_dir,\n        numbered_filename\n    )\n\n    plt.savefig(\n        output_path,\n        dpi=150,\n        bbox_inches=\"tight\"\n    )\n\n    plt.show()\n\n    print(\n        f\"{image_counter:02d}. {title}\"\n    )\n\n    print(\n        f\"   Saved: {output_path}\\n\"\n    )\n\n    plt.close()\n\n\n# ============================================================\n# 3. PREPARE COMMON DATA\n# ============================================================\n\ndaily_features[\"date\"] = pd.to_datetime(\n    daily_features[\"date\"]\n)\n\nforecast_df[\"date\"] = pd.to_datetime(\n    forecast_df[\"date\"]\n)\n\n# Make sure required series exist\nts = (\n    forecast_df\n    .set_index(\"date\")[\"total_sales\"]\n    .asfreq(\"D\")\n)\n\ndiff1 = ts.diff(1).dropna()\n\ndiff7 = ts.diff(7).dropna()\n\n\n# ============================================================\n# 4. IMAGE 01\n# DAILY SALES TREND\n# ============================================================\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    daily_sales[\"date\"],\n    daily_sales[\"total_sales\"]\n)\n\nplt.title(\n    \"Daily Total Sales Over Time\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.xticks(rotation=45)\n\nplt.tight_layout()\n\nsave_report_image(\n    \"Daily Total Sales Over Time\",\n    \"daily_sales_trend.png\"\n)\n\n\n# ============================================================\n# 5. IMAGE 02\n# DAILY SALES DISTRIBUTION\n# ============================================================\n\nplt.figure(figsize=(10, 5))\n\nplt.hist(\n    daily_sales[\"total_sales\"],\n    bins=50\n)\n\nplt.title(\n    \"Distribution of Daily Total Sales\"\n)\n\nplt.xlabel(\"Total Unit Sales\")\nplt.ylabel(\"Frequency\")\n\nplt.tight_layout()\n\nsave_report_image(\n    \"Distribution of Daily Total Sales\",\n    \"daily_sales_distribution.png\"\n)\n\n\n# ============================================================\n# 6. IMAGE 03\n# WEEKLY SALES TREND\n# ============================================================\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    weekly_sales[\"date\"],\n    weekly_sales[\"total_sales\"]\n)\n\nplt.title(\n    \"Weekly Total Sales Over Time\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.xticks(rotation=45)\n\nplt.tight_layout()\n\nsave_report_image(\n    \"Weekly Total Sales Over Time\",\n    \"weekly_sales_trend.png\"\n)\n\n\n# ============================================================\n# 7. IMAGE 04\n# MONTHLY SALES TREND\n# ============================================================\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    monthly_sales[\"date\"],\n    monthly_sales[\"total_sales\"]\n)\n\nplt.title(\n    \"Monthly Total Sales Over Time\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.xticks(rotation=45)\n\nplt.tight_layout()\n\nsave_report_image(\n    \"Monthly Total Sales Over Time\",\n    \"monthly_sales_trend.png\"\n)\n\n\n# ============================================================\n# 8. IMAGE 05\n# AVERAGE SALES BY MONTH\n# ============================================================\n\nplt.figure(figsize=(10, 5))\n\nplt.bar(\n    monthly_pattern.index,\n    monthly_pattern.values\n)\n\nplt.title(\n    \"Average Daily Sales by Month\"\n)\n\nplt.xlabel(\"Month\")\nplt.ylabel(\"Average Daily Sales\")\n\nplt.xticks(\n    range(1, 13)\n)\n\nplt.tight_layout()\n\nsave_report_image(\n    \"Average Daily Sales by Month\",\n    \"average_sales_by_month.png\"\n)\n\n\n# ============================================================\n# 9. IMAGE 06\n# AVERAGE SALES BY DAY OF WEEK\n# ============================================================\n\nplt.figure(figsize=(10, 5))\n\nplt.bar(\n    weekday_pattern[\"day_name\"],\n    weekday_pattern[\"total_sales\"]\n)\n\nplt.title(\n    \"Average Daily Sales by Day of Week\"\n)\n\nplt.xlabel(\"Day of Week\")\nplt.ylabel(\"Average Daily Sales\")\n\nplt.xticks(\n    rotation=45\n)\n\nplt.tight_layout()\n\nsave_report_image(\n    \"Average Daily Sales by Day of Week\",\n    \"average_sales_by_day_of_week.png\"\n)\n\n\n# ============================================================\n# 10. IMAGE 07\n# YEAR × MONTH HEATMAP\n# ============================================================\n\nplt.figure(figsize=(14, 6))\n\nsns.heatmap(\n    year_month_sales,\n    annot=True,\n    fmt=\".0f\"\n)\n\nplt.title(\n    \"Average Daily Sales by Year and Month\"\n)\n\nplt.xlabel(\"Month\")\nplt.ylabel(\"Year\")\n\nplt.tight_layout()\n\nsave_report_image(\n    \"Year by Month Sales Heatmap\",\n    \"year_month_sales_heatmap.png\"\n)\n\n\n# ============================================================\n# 11. IMAGE 08\n# HOLIDAY VS NON-HOLIDAY\n# ============================================================\n\nholiday_plot = (\n    daily_features\n    .dropna(subset=[\"total_sales\"])\n    .groupby(\"is_holiday\")[\"total_sales\"]\n    .mean()\n)\n\nholiday_labels = [\n    \"Non-Holiday\",\n    \"Holiday\"\n]\n\nholiday_values = [\n    holiday_plot.get(0, 0),\n    holiday_plot.get(1, 0)\n]\n\nplt.figure(figsize=(8, 5))\n\nplt.bar(\n    holiday_labels,\n    holiday_values\n)\n\nplt.title(\n    \"Average Sales: Holiday vs Non-Holiday\"\n)\n\nplt.xlabel(\"Day Type\")\nplt.ylabel(\"Average Daily Sales\")\n\nplt.tight_layout()\n\nsave_report_image(\n    \"Holiday vs Non-Holiday Sales\",\n    \"holiday_vs_nonholiday_sales.png\"\n)\n\n\n# ============================================================\n# 12. IMAGE 09\n# DAILY SALES + ROLLING TREND\n# ============================================================\n\ndaily_plot = daily_features.copy()\n\ndaily_plot[\"rolling_30d\"] = (\n    daily_plot[\"total_sales\"]\n    .rolling(\n        window=30,\n        min_periods=15\n    )\n    .mean()\n)\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    daily_plot[\"date\"],\n    daily_plot[\"total_sales\"],\n    alpha=0.35,\n    label=\"Daily Sales\"\n)\n\nplt.plot(\n    daily_plot[\"date\"],\n    daily_plot[\"rolling_30d\"],\n    linewidth=2,\n    label=\"30-Day Rolling Average\"\n)\n\nplt.title(\n    \"Daily Sales with 30-Day Rolling Trend\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\n\nplt.xticks(rotation=45)\n\nplt.tight_layout()\n\nsave_report_image(\n    \"Daily Sales with 30-Day Rolling Trend\",\n    \"daily_sales_30day_rolling.png\"\n)\n\n\n# ============================================================\n# 13. IMAGE 10\n# STATIONARITY — ORIGINAL\n# ============================================================\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    ts.index,\n    ts.values\n)\n\nplt.title(\n    \"Original Daily Sales Series\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.xticks(rotation=45)\n\nplt.tight_layout()\n\nsave_report_image(\n    \"Original Daily Sales Series\",\n    \"stationarity_original_series.png\"\n)\n\n\n# ============================================================\n# 14. IMAGE 11\n# FIRST DIFFERENCE\n# ============================================================\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    diff1.index,\n    diff1.values\n)\n\nplt.axhline(\n    0,\n    linestyle=\"--\"\n)\n\nplt.title(\n    \"First-Differenced Daily Sales\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Change in Sales\")\n\nplt.xticks(rotation=45)\n\nplt.tight_layout()\n\nsave_report_image(\n    \"First-Differenced Daily Sales\",\n    \"stationarity_first_difference.png\"\n)\n\n\n# ============================================================\n# 15. IMAGE 12\n# 7-DAY DIFFERENCE\n# ============================================================\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    diff7.index,\n    diff7.values\n)\n\nplt.axhline(\n    0,\n    linestyle=\"--\"\n)\n\nplt.title(\n    \"7-Day Seasonal Difference\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"7-Day Change in Sales\")\n\nplt.xticks(rotation=45)\n\nplt.tight_layout()\n\nsave_report_image(\n    \"7-Day Seasonal Difference\",\n    \"stationarity_7day_difference.png\"\n)\n\n\n# ============================================================\n# 16. IMAGE 13\n# ACF\n# ============================================================\n\nfig, ax = plt.subplots(\n    figsize=(14, 5)\n)\n\nplot_acf(\n    diff1,\n    lags=40,\n    ax=ax\n)\n\nax.set_title(\n    \"ACF — First-Differenced Daily Sales\"\n)\n\nplt.tight_layout()\n\nsave_report_image(\n    \"ACF Plot\",\n    \"acf_first_difference.png\"\n)\n\n\n# ============================================================\n# 17. IMAGE 14\n# PACF\n# ============================================================\n\nfig, ax = plt.subplots(\n    figsize=(14, 5)\n)\n\nplot_pacf(\n    diff1,\n    lags=40,\n    method=\"ywm\",\n    ax=ax\n)\n\nax.set_title(\n    \"PACF — First-Differenced Daily Sales\"\n)\n\nplt.tight_layout()\n\nsave_report_image(\n    \"PACF Plot\",\n    \"pacf_first_difference.png\"\n)\n\n\n# ============================================================\n# 18. IMAGE 15\n# SEASONAL ACF\n# ============================================================\n\nfig, ax = plt.subplots(\n    figsize=(14, 5)\n)\n\nplot_acf(\n    diff7,\n    lags=60,\n    ax=ax\n)\n\nax.set_title(\n    \"ACF — 7-Day Seasonal Difference\"\n)\n\nplt.tight_layout()\n\nsave_report_image(\n    \"Seasonal ACF Plot\",\n    \"acf_7day_difference.png\"\n)\n\n\n# ============================================================\n# 19. IMAGE 16\n# SEASONAL PACF\n# ============================================================\n\nfig, ax = plt.subplots(\n    figsize=(14, 5)\n)\n\nplot_pacf(\n    diff7,\n    lags=60,\n    method=\"ywm\",\n    ax=ax\n)\n\nax.set_title(\n    \"PACF — 7-Day Seasonal Difference\"\n)\n\nplt.tight_layout()\n\nsave_report_image(\n    \"Seasonal PACF Plot\",\n    \"pacf_7day_difference.png\"\n)\n\n\n# ============================================================\n# 20. IMAGE 17\n# VALIDATION — SEASONAL NAIVE\n# ============================================================\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    validation_ts.index,\n    validation_ts.values,\n    label=\"Actual\"\n)\n\nplt.plot(\n    validation_ts.index,\n    seasonal_naive_predictions.values,\n    label=\"Seasonal Naive\"\n)\n\nplt.title(\n    \"Actual vs Seasonal Naive — Validation\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\n\nplt.xticks(rotation=45)\n\nplt.tight_layout()\n\nsave_report_image(\n    \"Actual vs Seasonal Naive — Validation\",\n    \"validation_actual_vs_seasonal_naive.png\"\n)\n\n\n# ============================================================\n# 21. IMAGE 18\n# VALIDATION — ARIMA\n# ============================================================\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    validation_ts.index,\n    validation_ts.values,\n    label=\"Actual\"\n)\n\nplt.plot(\n    validation_ts.index,\n    arima_predictions.values,\n    label=\"ARIMA\"\n)\n\nplt.title(\n    \"Actual vs ARIMA — Validation\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\n\nplt.xticks(rotation=45)\n\nplt.tight_layout()\n\nsave_report_image(\n    \"Actual vs ARIMA — Validation\",\n    \"validation_actual_vs_arima.png\"\n)\n\n\n# ============================================================\n# 22. IMAGE 19\n# VALIDATION — SARIMA\n# ============================================================\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    validation_ts.index,\n    validation_ts.values,\n    label=\"Actual\"\n)\n\nplt.plot(\n    validation_ts.index,\n    sarima_predictions.values,\n    label=\"SARIMA\"\n)\n\nplt.title(\n    \"Actual vs SARIMA — Validation\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\n\nplt.xticks(rotation=45)\n\nplt.tight_layout()\n\nsave_report_image(\n    \"Actual vs SARIMA — Validation\",\n    \"validation_actual_vs_sarima.png\"\n)\n\n\n# ============================================================\n# 23. IMAGE 20\n# VALIDATION — EXPONENTIAL SMOOTHING\n# ============================================================\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    validation_ts.index,\n    validation_ts.values,\n    label=\"Actual\"\n)\n\nplt.plot(\n    validation_ts.index,\n    exp_predictions.values,\n    label=\"Exponential Smoothing\"\n)\n\nplt.title(\n    \"Actual vs Exponential Smoothing — Validation\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\n\nplt.xticks(rotation=45)\n\nplt.tight_layout()\n\nsave_report_image(\n    \"Actual vs Exponential Smoothing — Validation\",\n    \"validation_actual_vs_exponential_smoothing.png\"\n)\n\n\n# ============================================================\n# 24. IMAGE 21\n# VALIDATION — PROPHET\n# ============================================================\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    validation_ts.index,\n    validation_ts.values,\n    label=\"Actual\"\n)\n\nplt.plot(\n    validation_ts.index,\n    prophet_predictions.values,\n    label=\"Prophet\"\n)\n\nplt.title(\n    \"Actual vs Prophet — Validation\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\n\nplt.xticks(rotation=45)\n\nplt.tight_layout()\n\nsave_report_image(\n    \"Actual vs Prophet — Validation\",\n    \"validation_actual_vs_prophet.png\"\n)\n\n\n# ============================================================\n# 25. IMAGE 22\n# VALIDATION — LSTM\n# ============================================================\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    validation_ts.index,\n    validation_ts.values,\n    label=\"Actual\"\n)\n\nplt.plot(\n    validation_ts.index,\n    lstm_predictions.values,\n    label=\"LSTM\"\n)\n\nplt.title(\n    \"Actual vs LSTM — Validation\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\n\nplt.xticks(rotation=45)\n\nplt.tight_layout()\n\nsave_report_image(\n    \"Actual vs LSTM — Validation\",\n    \"validation_actual_vs_lstm.png\"\n)\n\n\n# ============================================================\n# 26. IMAGE 23\n# VALIDATION — MAIN MODELS TOGETHER\n# ============================================================\n\nplt.figure(figsize=(18, 7))\n\nplt.plot(\n    validation_ts.index,\n    validation_ts.values,\n    label=\"Actual\"\n)\n\nplt.plot(\n    validation_ts.index,\n    sarima_predictions.values,\n    label=\"SARIMA\"\n)\n\nplt.plot(\n    validation_ts.index,\n    exp_predictions.values,\n    label=\"Exponential Smoothing\"\n)\n\nplt.plot(\n    validation_ts.index,\n    prophet_predictions.values,\n    label=\"Prophet\"\n)\n\nplt.plot(\n    validation_ts.index,\n    lstm_predictions.values,\n    label=\"LSTM\"\n)\n\nplt.title(\n    \"Actual vs Major Forecasting Models\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\n\nplt.xticks(rotation=45)\n\nplt.tight_layout()\n\nsave_report_image(\n    \"Actual vs SARIMA vs Exponential Smoothing vs Prophet vs LSTM\",\n    \"validation_all_main_models.png\"\n)\n\n\n# ============================================================\n# 27. IMAGE 24\n# BACKTEST MAE BY FOLD\n# ============================================================\n\nsarima_fold_results = (\n    backtest_results_df[\n        backtest_results_df[\"Model Type\"] == \"SARIMA\"\n    ]\n)\n\nseasonal_fold_results = (\n    backtest_results_df[\n        backtest_results_df[\"Model Type\"] == \"Seasonal Naive\"\n    ]\n)\n\nplt.figure(figsize=(10, 5))\n\nplt.plot(\n    sarima_fold_results[\"Fold\"],\n    sarima_fold_results[\"MAE\"],\n    marker=\"o\",\n    label=\"SARIMA\"\n)\n\nplt.plot(\n    seasonal_fold_results[\"Fold\"],\n    seasonal_fold_results[\"MAE\"],\n    marker=\"o\",\n    label=\"Seasonal Naive\"\n)\n\nplt.title(\n    \"Backtesting MAE by Fold\"\n)\n\nplt.xlabel(\"Backtest Fold\")\nplt.ylabel(\"MAE\")\n\nplt.legend()\n\nplt.tight_layout()\n\nsave_report_image(\n    \"Backtesting MAE by Fold\",\n    \"backtesting_mae_by_fold.png\"\n)\n\n\n# ============================================================\n# 28. IMAGE 25\n# BACKTEST FOLD 1\n# ============================================================\n\nfirst_fold = (\n    backtest_predictions_df[\n        backtest_predictions_df[\"fold\"] == 1\n    ]\n)\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    first_fold[\"date\"],\n    first_fold[\"actual\"],\n    label=\"Actual\"\n)\n\nplt.plot(\n    first_fold[\"date\"],\n    first_fold[\"sarima_prediction\"],\n    label=\"SARIMA\"\n)\n\nplt.plot(\n    first_fold[\"date\"],\n    first_fold[\"seasonal_naive_prediction\"],\n    label=\"Seasonal Naive\"\n)\n\nplt.title(\n    \"Backtest Fold 1 — Actual vs Forecast\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\n\nplt.xticks(rotation=45)\n\nplt.tight_layout()\n\nsave_report_image(\n    \"Backtest Fold 1 Actual vs Forecast\",\n    \"backtest_fold1_actual_vs_forecast.png\"\n)\n\n\n# ============================================================\n# 29. IMAGE 26\n# FINAL TEST — TUNED SARIMA\n# ============================================================\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    test_ts.index,\n    test_ts.values,\n    label=\"Actual\"\n)\n\nplt.plot(\n    test_predictions.index,\n    test_predictions.values,\n    label=\"Tuned SARIMA\"\n)\n\nplt.title(\n    \"Final Test — Actual vs Tuned SARIMA\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\n\nplt.xticks(rotation=45)\n\nplt.tight_layout()\n\nsave_report_image(\n    \"Final Test Actual vs Tuned SARIMA\",\n    \"final_test_actual_vs_sarima.png\"\n)\n\n\n# ============================================================\n# 30. IMAGE 27\n# FINAL TEST — SEASONAL NAIVE\n# ============================================================\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    test_ts.index,\n    test_ts.values,\n    label=\"Actual\"\n)\n\nplt.plot(\n    seasonal_naive_test_predictions.index,\n    seasonal_naive_test_predictions.values,\n    label=\"Seasonal Naive\"\n)\n\nplt.title(\n    \"Final Test — Actual vs Seasonal Naive\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\n\nplt.xticks(rotation=45)\n\nplt.tight_layout()\n\nsave_report_image(\n    \"Final Test Actual vs Seasonal Naive\",\n    \"final_test_actual_vs_seasonal_naive.png\"\n)\n\n\n# ============================================================\n# 31. IMAGE 28\n# FINAL TEST — BOTH MODELS\n# ============================================================\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    test_ts.index,\n    test_ts.values,\n    label=\"Actual\"\n)\n\nplt.plot(\n    test_predictions.index,\n    test_predictions.values,\n    label=\"Tuned SARIMA\"\n)\n\nplt.plot(\n    seasonal_naive_test_predictions.index,\n    seasonal_naive_test_predictions.values,\n    label=\"Seasonal Naive\"\n)\n\nplt.title(\n    \"Final Test — Actual vs Tuned SARIMA vs Seasonal Naive\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\n\nplt.xticks(rotation=45)\n\nplt.tight_layout()\n\nsave_report_image(\n    \"Final Test Actual vs Tuned SARIMA vs Seasonal Naive\",\n    \"final_test_all_models.png\"\n)\n\n\n# ============================================================\n# 32. IMAGE 29\n# MODEL MAE COMPARISON\n# ============================================================\n\nplt.figure(figsize=(12, 6))\n\nplt.bar(\n    final_model_comparison[\"Model\"],\n    final_model_comparison[\"MAE\"]\n)\n\nplt.title(\n    \"Forecasting Model Comparison — MAE\"\n)\n\nplt.xlabel(\"Model\")\nplt.ylabel(\"MAE\")\n\nplt.xticks(\n    rotation=45,\n    ha=\"right\"\n)\n\nplt.tight_layout()\n\nsave_report_image(\n    \"Forecasting Model Comparison — MAE\",\n    \"model_comparison_mae.png\"\n)\n\n\n# ============================================================\n# 33. IMAGE 30\n# FINAL 6-MONTH FORECAST\n# ============================================================\n\nrecent_history = final_series.tail(180)\n\nplt.figure(figsize=(17, 7))\n\nplt.plot(\n    recent_history.index,\n    recent_history.values,\n    label=\"Historical Sales\"\n)\n\nplt.plot(\n    forecast_6m.index,\n    forecast_6m[\"forecast\"],\n    label=\"6-Month Forecast\"\n)\n\nplt.fill_between(\n    forecast_6m.index,\n    forecast_6m[\"lower_95_inventory\"],\n    forecast_6m[\"upper_95_inventory\"],\n    alpha=0.2,\n    label=\"95% Confidence Interval\"\n)\n\nplt.axvline(\n    final_series.index.max(),\n    linestyle=\"--\",\n    label=\"Forecast Start\"\n)\n\nplt.title(\n    \"6-Month Demand Forecast with 95% Confidence Interval\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\n\nplt.xticks(rotation=45)\n\nplt.tight_layout()\n\nsave_report_image(\n    \"6-Month Forecast with 95% Confidence Interval\",\n    \"forecast_6_month_confidence_interval.png\"\n)\n\n\n# ============================================================\n# 34. IMAGE 31\n# FINAL 12-MONTH FORECAST\n# ============================================================\n\nplt.figure(figsize=(18, 7))\n\nplt.plot(\n    recent_history.index,\n    recent_history.values,\n    label=\"Historical Sales\"\n)\n\nplt.plot(\n    forecast_output.index,\n    forecast_output[\"forecast\"],\n    label=\"12-Month Forecast\"\n)\n\nplt.fill_between(\n    forecast_output.index,\n    forecast_output[\"lower_95_inventory\"],\n    forecast_output[\"upper_95_inventory\"],\n    alpha=0.2,\n    label=\"95% Confidence Interval\"\n)\n\nplt.axvline(\n    final_series.index.max(),\n    linestyle=\"--\",\n    label=\"Forecast Start\"\n)\n\nplt.title(\n    \"Historical Sales and 12-Month SARIMA Forecast\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\n\nplt.xticks(rotation=45)\n\nplt.tight_layout()\n\nsave_report_image(\n    \"12-Month Forecast with 95% Confidence Interval\",\n    \"forecast_12_month_confidence_interval.png\"\n)\n\n\n# ============================================================\n# 35. IMAGE 32\n# MONTHLY FORECAST\n# ============================================================\n\nplt.figure(figsize=(14, 6))\n\nplt.bar(\n    monthly_forecast[\"date\"].astype(str),\n    monthly_forecast[\"forecast_units\"]\n)\n\nplt.title(\n    \"Forecast Monthly Demand\"\n)\n\nplt.xlabel(\"Month\")\nplt.ylabel(\"Forecast Unit Sales\")\n\nplt.xticks(\n    rotation=60\n)\n\nplt.tight_layout()\n\nsave_report_image(\n    \"Forecast Monthly Demand\",\n    \"monthly_forecast_demand.png\"\n)\n\n\n# ============================================================\n# 36. IMAGE 33\n# FORECAST UNCERTAINTY\n# ============================================================\n\nforecast_interval_width = (\n    forecast_output[\"upper_95_inventory\"]\n    -\n    forecast_output[\"lower_95_inventory\"]\n)\n\nplt.figure(figsize=(16, 6))\n\nplt.plot(\n    forecast_output.index,\n    forecast_interval_width\n)\n\nplt.title(\n    \"Forecast Uncertainty Over 12 Months\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"95% Confidence Interval Width\")\n\nplt.xticks(rotation=45)\n\nplt.tight_layout()\n\nsave_report_image(\n    \"Forecast Uncertainty Over 12 Months\",\n    \"forecast_uncertainty.png\"\n)\n\n\n# ============================================================\n# 37. IMAGE 34\n# INVENTORY PLANNING UNCERTAINTY\n# ============================================================\n\nplt.figure(figsize=(17, 7))\n\nplt.plot(\n    inventory_df.index,\n    inventory_df[\"forecast\"],\n    label=\"Expected Demand\"\n)\n\nplt.fill_between(\n    inventory_df.index,\n    inventory_df[\"lower_95\"],\n    inventory_df[\"upper_95\"],\n    alpha=0.2,\n    label=\"95% Forecast Interval\"\n)\n\nplt.title(\n    \"Forecast Demand and Uncertainty for Inventory Planning\"\n)\n\nplt.xlabel(\"Date\")\nplt.ylabel(\"Total Unit Sales\")\n\nplt.legend()\n\nplt.xticks(rotation=45)\n\nplt.tight_layout()\n\nsave_report_image(\n    \"Forecast Demand and Uncertainty for Inventory Planning\",\n    \"inventory_demand_uncertainty.png\"\n)\n\n\n# ============================================================\n# 38. IMAGE 35\n# PRODUCT-LEVEL MODEL ACCURACY\n# ============================================================\n\nproduct_plot_data = (\n    product_accuracy[\n        [\n            \"item_nbr\",\n            \"Model\",\n            \"MAE\"\n        ]\n    ]\n    .copy()\n)\n\nproduct_plot_data[\"Model Type\"] = (\n    product_plot_data[\"Model\"]\n    .str.replace(\n        r\" - \\d+\",\n        \"\",\n        regex=True\n    )\n)\n\nproduct_mae_pivot = (\n    product_plot_data\n    .pivot(\n        index=\"item_nbr\",\n        columns=\"Model Type\",\n        values=\"MAE\"\n    )\n)\n\nplt.figure(figsize=(12, 6))\n\nproduct_mae_pivot.plot(\n    kind=\"bar\",\n    ax=plt.gca()\n)\n\nplt.title(\n    \"Product-Level Model Accuracy — MAE\"\n)\n\nplt.xlabel(\"Product\")\nplt.ylabel(\"MAE\")\n\nplt.xticks(\n    rotation=45\n)\n\nplt.tight_layout()\n\nsave_report_image(\n    \"Product-Level Model Accuracy Comparison\",\n    \"product_level_model_accuracy.png\"\n)\n\n\n# ============================================================\n# 39. IMAGES 36–40\n# PRODUCT-LEVEL ACTUAL VS PREDICTED\n# ============================================================\n\nfor product_number, item_id in enumerate(\n    top_product_ids,\n    start=1\n):\n\n    if item_id not in product_predictions:\n        continue\n\n    result = product_predictions[item_id]\n\n    try:\n\n        family = (\n            top_5_products\n            .loc[\n                top_5_products[\"item_nbr\"] == item_id,\n                \"family\"\n            ]\n            .iloc[0]\n        )\n\n    except Exception:\n\n        family = \"Unknown\"\n\n\n    plt.figure(figsize=(16, 6))\n\n    plt.plot(\n        result[\"actual\"].index,\n        result[\"actual\"].values,\n        label=\"Actual\"\n    )\n\n    plt.plot(\n        result[\"sarima\"].index,\n        result[\"sarima\"].values,\n        label=\"SARIMA\"\n    )\n\n    plt.plot(\n        result[\"seasonal_naive\"].index,\n        result[\"seasonal_naive\"].values,\n        label=\"Seasonal Naive\"\n    )\n\n    plt.title(\n        f\"Product {item_id} ({family}) — \"\n        \"Actual vs Predicted\"\n    )\n\n    plt.xlabel(\"Date\")\n    plt.ylabel(\"Unit Sales\")\n\n    plt.legend()\n\n    plt.xticks(rotation=45)\n\n    plt.tight_layout()\n\n    save_report_image(\n        f\"Product {product_number} — \"\n        f\"{item_id} Actual vs Predicted\",\n        f\"product_{product_number}_{item_id}_actual_vs_predicted.png\"\n    )\n\n\n# ============================================================\n# 40. FINAL IMAGE LIST\n# ============================================================\n\nprint(\"=\" * 80)\nprint(\"ALL REPORT IMAGES GENERATED\")\nprint(\"=\" * 80)\n\nprint(\n    f\"\\nTotal images generated: \"\n    f\"{image_counter}\"\n)\n\nprint(\n    f\"\\nImage folder:\\n{image_dir}\"\n)\n\nprint(\"\\nOrdered image files:\")\n\nfor filename in sorted(\n    os.listdir(image_dir)\n):\n\n    print(\n        os.path.join(\n            image_dir,\n            filename\n        )\n    )","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:45:34.135043Z","iopub.execute_input":"2026-09-17T10:45:34.136096Z","iopub.status.idle":"2026-09-17T10:46:00.540337Z","shell.execute_reply.started":"2026-09-17T10:45:34.136058Z","shell.execute_reply":"2026-09-17T10:46:00.539226Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# STEP 35 — CREATE ONE SELF-CONTAINED REPORT FILE\n# ============================================================\n#\n# Output:\n# /kaggle/working/Favorita_Demand_Forecasting_Final_Report.html\n#\n# This HTML contains:\n#   - Full Markdown report\n#   - All generated plots embedded inside the file\n#   - No dependency on separate image files\n#\n# You can:\n#   1. Open the HTML in a browser\n#   2. Review the complete report\n#   3. Print / Save as PDF\n#\n# ============================================================\n\nfrom pathlib import Path\nimport base64\nimport re\nimport html\n\n\n# ============================================================\n# 1. PATHS\n# ============================================================\n\nreport_md_path = Path(\n    \"/kaggle/working/Favorita_Demand_Forecasting_Report.md\"\n)\n\nimage_dir = Path(\n    \"/kaggle/working/final_report_images\"\n)\n\noutput_html_path = Path(\n    \"/kaggle/working/\"\n    \"Favorita_Demand_Forecasting_Final_Report.html\"\n)\n\n\n# ============================================================\n# 2. CHECK FILES\n# ============================================================\n\nif not report_md_path.exists():\n    raise FileNotFoundError(\n        f\"Report not found:\\n{report_md_path}\"\n    )\n\nif not image_dir.exists():\n    raise FileNotFoundError(\n        f\"Image directory not found:\\n{image_dir}\"\n    )\n\n\n# ============================================================\n# 3. READ MARKDOWN REPORT\n# ============================================================\n\nmarkdown_text = report_md_path.read_text(\n    encoding=\"utf-8\"\n)\n\n\n# ============================================================\n# 4. IMAGE FILE LIST\n# ============================================================\n\nimage_files = sorted(\n    image_dir.glob(\"*.png\")\n)\n\nprint(\"=\" * 80)\nprint(\"REPORT + IMAGE COMBINER\")\nprint(\"=\" * 80)\n\nprint(\n    f\"\\nMarkdown report:\\n{report_md_path}\"\n)\n\nprint(\n    f\"\\nImages found: {len(image_files)}\"\n)\n\nfor image in image_files:\n    print(\n        \" -\",\n        image.name\n    )\n\n\n# ============================================================\n# 5. MAP REPORT PLACEHOLDERS TO IMAGE FILES\n# ============================================================\n\nimage_map = {\n\n    # --------------------------------------------------------\n    # EDA\n    # --------------------------------------------------------\n\n    \"Daily Total Sales Over Time\":\n        \"01_daily_sales_trend.png\",\n\n    \"Distribution of Daily Total Sales\":\n        \"02_daily_sales_distribution.png\",\n\n    \"Weekly Total Sales Over Time\":\n        \"03_weekly_sales_trend.png\",\n\n    \"Monthly Total Sales Over Time\":\n        \"04_monthly_sales_trend.png\",\n\n    \"Average Daily Sales by Month\":\n        \"05_average_sales_by_month.png\",\n\n    \"Average Daily Sales by Day of Week\":\n        \"06_average_sales_by_day_of_week.png\",\n\n    \"Year by Month Sales Heatmap\":\n        \"07_year_month_sales_heatmap.png\",\n\n    \"Holiday vs Non-Holiday Sales\":\n        \"08_holiday_vs_nonholiday_sales.png\",\n\n    # --------------------------------------------------------\n    # Stationarity\n    # --------------------------------------------------------\n\n    \"Original Daily Sales Series\":\n        \"10_stationarity_original_series.png\",\n\n    \"First-Differenced Daily Sales\":\n        \"11_stationarity_first_difference.png\",\n\n    \"7-Day Seasonal Difference\":\n        \"12_stationarity_7day_difference.png\",\n\n    # --------------------------------------------------------\n    # ACF / PACF\n    # --------------------------------------------------------\n\n    \"ACF Plot\":\n        \"13_acf_first_difference.png\",\n\n    \"PACF Plot\":\n        \"14_pacf_first_difference.png\",\n\n    \"Seasonal ACF Plot\":\n        \"15_acf_7day_difference.png\",\n\n    \"Seasonal PACF Plot\":\n        \"16_pacf_7day_difference.png\",\n\n    # --------------------------------------------------------\n    # Validation\n    # --------------------------------------------------------\n\n    \"Actual vs Seasonal Naive — Validation\":\n        \"17_validation_actual_vs_seasonal_naive.png\",\n\n    \"Actual vs ARIMA — Validation\":\n        \"18_validation_actual_vs_arima.png\",\n\n    \"Actual vs SARIMA — Validation\":\n        \"19_validation_actual_vs_sarima.png\",\n\n    \"Actual vs Exponential Smoothing — Validation\":\n        \"20_validation_actual_vs_exponential_smoothing.png\",\n\n    \"Actual vs Prophet — Validation\":\n        \"21_validation_actual_vs_prophet.png\",\n\n    \"Actual vs LSTM — Validation\":\n        \"22_validation_actual_vs_lstm.png\",\n\n    \"Actual vs SARIMA vs Exponential Smoothing vs Prophet vs LSTM\":\n        \"23_validation_all_main_models.png\",\n\n    # --------------------------------------------------------\n    # Backtesting\n    # --------------------------------------------------------\n\n    \"Backtesting MAE by Fold\":\n        \"24_backtesting_mae_by_fold.png\",\n\n    \"Backtest Fold 1 Actual vs Forecast\":\n        \"25_backtest_fold1_actual_vs_forecast.png\",\n\n    # --------------------------------------------------------\n    # Final test\n    # --------------------------------------------------------\n\n    \"Final Test Actual vs Tuned SARIMA\":\n        \"26_final_test_actual_vs_sarima.png\",\n\n    \"Final Test Actual vs Seasonal Naive\":\n        \"27_final_test_actual_vs_seasonal_naive.png\",\n\n    \"Final Test Actual vs Tuned SARIMA vs Seasonal Naive\":\n        \"28_final_test_all_models.png\",\n\n    # --------------------------------------------------------\n    # Forecasting\n    # --------------------------------------------------------\n\n    \"6-Month Forecast with 95% Confidence Interval\":\n        \"30_forecast_6_month_confidence_interval.png\",\n\n    \"12-Month Forecast with 95% Confidence Interval\":\n        \"31_forecast_12_month_confidence_interval.png\",\n\n    \"Forecast Monthly Demand\":\n        \"32_monthly_forecast_demand.png\",\n\n    \"Forecast Uncertainty Over 12 Months\":\n        \"33_forecast_uncertainty.png\",\n\n    \"Forecast Demand and Uncertainty for Inventory Planning\":\n        \"34_inventory_demand_uncertainty.png\",\n\n    # --------------------------------------------------------\n    # Product level\n    # --------------------------------------------------------\n\n    \"Product-Level Model Accuracy Comparison\":\n        \"35_product_level_model_accuracy.png\",\n}\n\n\n# ============================================================\n# 6. PRODUCT IMAGE MAPPING\n# ============================================================\n\nproduct_images = sorted(\n    image_dir.glob(\n        \"3[6-9]_product_*_actual_vs_predicted.png\"\n    )\n)\n\n# Also capture image 40\nproduct_images += sorted(\n    image_dir.glob(\n        \"40_product_*_actual_vs_predicted.png\"\n    )\n)\n\n# More robust fallback:\nif len(product_images) < 5:\n\n    all_product_images = sorted(\n        image_dir.glob(\n            \"*_product_*_actual_vs_predicted.png\"\n        )\n    )\n\n    product_images = all_product_images\n\n\nprint(\n    f\"\\nProduct-level images found: \"\n    f\"{len(product_images)}\"\n)\n\n\n# ============================================================\n# 7. CONVERT IMAGE TO BASE64 HTML\n# ============================================================\n\ndef image_to_html(\n    image_path,\n    title\n):\n\n    image_bytes = image_path.read_bytes()\n\n    encoded = base64.b64encode(\n        image_bytes\n    ).decode(\"utf-8\")\n\n    return f\"\"\"\n<figure class=\"report-figure\">\n    <img\n        src=\"data:image/png;base64,{encoded}\"\n        alt=\"{html.escape(title)}\"\n    />\n    <figcaption>\n        {html.escape(title)}\n    </figcaption>\n</figure>\n\"\"\"\n\n\n# ============================================================\n# 8. REPLACE STANDARD IMAGE PLACEHOLDERS\n# ============================================================\n\nembedded_count = 0\n\nfor title, filename in image_map.items():\n\n    placeholder = (\n        f\"[image - {title}]\"\n    )\n\n    image_path = (\n        image_dir / filename\n    )\n\n    if image_path.exists():\n\n        replacement = image_to_html(\n            image_path,\n            title\n        )\n\n        markdown_text = (\n            markdown_text.replace(\n                placeholder,\n                replacement\n            )\n        )\n\n        embedded_count += 1\n\n    else:\n\n        print(\n            f\"WARNING: Image missing for: \"\n            f\"{title}\"\n        )\n\n\n# ============================================================\n# 9. REPLACE PRODUCT PLACEHOLDERS\n# ============================================================\n\nfor i in range(1, 6):\n\n    placeholder = (\n        f\"[image - Product {i} Actual vs Predicted]\"\n    )\n\n    if i <= len(product_images):\n\n        product_image = product_images[i - 1]\n\n        replacement = image_to_html(\n            product_image,\n            f\"Product {i} — Actual vs Predicted\"\n        )\n\n        markdown_text = (\n            markdown_text.replace(\n                placeholder,\n                replacement\n            )\n        )\n\n        embedded_count += 1\n\n    else:\n\n        print(\n            f\"WARNING: Product {i} image not found.\"\n        )\n\n\n# ============================================================\n# 10. HANDLE THE ROLLING TREND IMAGE\n# ============================================================\n\nrolling_image = (\n    image_dir /\n    \"09_daily_sales_30day_rolling.png\"\n)\n\nrolling_placeholder = (\n    \"[image - Daily Sales with 30-Day Rolling Trend]\"\n)\n\nif rolling_image.exists():\n\n    # Add this visualization immediately after\n    # the EDA visualization section.\n    rolling_html = image_to_html(\n        rolling_image,\n        \"Daily Sales with 30-Day Rolling Trend\"\n    )\n\n    # If the placeholder isn't in the report,\n    # we'll append it later in the appendix.\n    if rolling_placeholder in markdown_text:\n\n        markdown_text = (\n            markdown_text.replace(\n                rolling_placeholder,\n                rolling_html\n            )\n        )\n\n\n# ============================================================\n# 11. OPTIONAL: CREATE APPENDIX FOR ALL EXTRA IMAGES\n# ============================================================\n#\n# Some generated images may not have an explicit placeholder\n# in the Markdown report. We don't want to lose them.\n#\n# These include:\n#   - Rolling trend\n#   - Stationarity plots if not used\n#   - Model comparison chart\n#\n# ============================================================\n\nused_files = set(\n    image_map.values()\n)\n\nfor image_path in product_images:\n    used_files.add(\n        image_path.name\n    )\n\nextra_images = [\n    image\n    for image in image_files\n    if image.name not in used_files\n]\n\n\nif extra_images:\n\n    appendix = \"\"\"\n\n# Appendix — Additional Generated Visualizations\n\nThe following visualizations were generated during the analysis\nand are included here for completeness.\n\n\"\"\"\n\n    for image_path in extra_images:\n\n        title = (\n            image_path\n            .stem\n            .replace(\"_\", \" \")\n            .title()\n        )\n\n        appendix += image_to_html(\n            image_path,\n            title\n        )\n\n    markdown_text += appendix\n\n\n# ============================================================\n# 12. CONVERT MARKDOWN TO HTML\n# ============================================================\n\ntry:\n\n    import markdown\n\n    html_body = markdown.markdown(\n        markdown_text,\n        extensions=[\n            \"tables\",\n            \"fenced_code\",\n            \"toc\"\n        ]\n    )\n\nexcept ImportError:\n\n    print(\n        \"\\nMarkdown package not found.\"\n    )\n\n    print(\n        \"Installing markdown...\"\n    )\n\n    !pip install markdown -q\n\n    import markdown\n\n    html_body = markdown.markdown(\n        markdown_text,\n        extensions=[\n            \"tables\",\n            \"fenced_code\",\n            \"toc\"\n        ]\n    )\n\n\n# ============================================================\n# 13. CREATE PROFESSIONAL HTML DOCUMENT\n# ============================================================\n\nhtml_document = f\"\"\"\n<!DOCTYPE html>\n\n<html lang=\"en\">\n\n<head>\n\n<meta charset=\"UTF-8\">\n\n<meta name=\"viewport\"\n      content=\"width=device-width,\n               initial-scale=1.0\">\n\n<title>\nFavorita Grocery Demand Forecasting\n</title>\n\n\n<style>\n\nbody {{\n    font-family:\n        Arial,\n        Helvetica,\n        sans-serif;\n\n    line-height: 1.6;\n\n    color: #222;\n\n    max-width: 1200px;\n\n    margin: 0 auto;\n\n    padding: 40px;\n}}\n\nh1 {{\n    font-size: 32px;\n\n    border-bottom:\n        3px solid #222;\n\n    padding-bottom: 12px;\n}}\n\nh2 {{\n    font-size: 25px;\n\n    margin-top: 45px;\n\n    border-bottom:\n        1px solid #ddd;\n\n    padding-bottom: 8px;\n}}\n\nh3 {{\n    font-size: 20px;\n\n    margin-top: 30px;\n}}\n\ntable {{\n    border-collapse: collapse;\n\n    width: 100%;\n\n    margin: 20px 0;\n}}\n\nth,\ntd {{\n    border: 1px solid #ddd;\n\n    padding: 9px;\n\n    text-align: left;\n}}\n\nth {{\n    background: #f2f2f2;\n\n    font-weight: bold;\n}}\n\ntr:nth-child(even) {{\n    background: #fafafa;\n}}\n\n.report-figure {{\n    margin: 30px auto;\n\n    text-align: center;\n\n    page-break-inside:\n        avoid;\n}}\n\n.report-figure img {{\n    max-width: 100%;\n\n    height: auto;\n\n    display: block;\n\n    margin: 0 auto;\n\n    border:\n        1px solid #ddd;\n}}\n\n.report-figure figcaption {{\n    margin-top: 8px;\n\n    font-weight: bold;\n\n    font-size: 14px;\n\n    color: #444;\n}}\n\ncode {{\n    background: #f4f4f4;\n\n    padding: 2px 5px;\n\n    border-radius: 3px;\n}}\n\npre {{\n    background: #f4f4f4;\n\n    padding: 15px;\n\n    overflow-x: auto;\n}}\n\nhr {{\n    border: 0;\n\n    border-top:\n        1px solid #ddd;\n\n    margin: 35px 0;\n}}\n\n@media print {{\n\n    body {{\n        max-width: none;\n\n        padding: 20px;\n    }}\n\n    h1,\n    h2,\n    h3 {{\n        page-break-after:\n            avoid;\n    }}\n\n    .report-figure {{\n        page-break-inside:\n            avoid;\n    }}\n\n    table {{\n        page-break-inside:\n            avoid;\n    }}\n\n}}\n\n</style>\n\n</head>\n\n\n<body>\n\n{html_body}\n\n</body>\n\n</html>\n\"\"\"\n\n\n# ============================================================\n# 14. SAVE SINGLE FILE\n# ============================================================\n\noutput_html_path.write_text(\n    html_document,\n    encoding=\"utf-8\"\n)\n\n\n# ============================================================\n# 15. FINAL CHECK\n# ============================================================\n\nfile_size_mb = (\n    output_html_path.stat().st_size\n    /\n    (1024 ** 2)\n)\n\n\nprint(\"\\n\" + \"=\" * 80)\nprint(\"SINGLE-FILE REPORT CREATED\")\nprint(\"=\" * 80)\n\nprint(\n    \"\\nFile:\"\n)\n\nprint(\n    output_html_path\n)\n\nprint(\n    f\"\\nFile size: \"\n    f\"{file_size_mb:.2f} MB\"\n)\n\nprint(\n    f\"\\nEmbedded standard images: \"\n    f\"{embedded_count}\"\n)\n\nprint(\n    f\"Total PNG images available: \"\n    f\"{len(image_files)}\"\n)\n\nprint(\n    \"\\nEverything is embedded inside this ONE HTML file.\"\n)\n\nprint(\n    \"\\nOpen the HTML file in your browser, \"\n    \"then use Print → Save as PDF \"\n    \"to create a single PDF version.\"\n)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-17T10:47:27.059089Z","iopub.execute_input":"2026-09-17T10:47:27.060158Z","iopub.status.idle":"2026-09-17T10:47:28.432377Z","shell.execute_reply.started":"2026-09-17T10:47:27.060116Z","shell.execute_reply":"2026-09-17T10:47:28.43129Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}