diff --git a/README.md b/README.md
index 2e58fea..b13755e 100644
--- a/README.md
+++ b/README.md
@@ -237,6 +237,32 @@ While both models demonstrated similar overall accuracy, the LSTM model's abilit
 
 For future work, combining the strengths of ARIMA and LSTM could be a promising approach. While ARIMA is effective for modeling linear trends and seasonality, LSTM excels in capturing complex patterns and sudden changes. A hybrid model leveraging these complementary strengths, along with the insights gained from sensitivity analysis, may yield even better forecasting results.
 
+## Recommendations for Real-World Application
+
+The PSO-LSTM model can be effectively used in industrial wastewater treatment to optimize processes and enhance decision-making. Here are some key recommendations:
+
+### 1. **Real-Time Monitoring and Control**
+   - **Integration with Monitoring Systems:** Use the model for real-time predictions of heavy metal concentrations to dynamically adjust treatment parameters, such as chemical dosage and flow rate.
+   - **Automated Control Systems:** Enable automated process adjustments to maintain compliance with environmental standards and optimize treatment efficiency.
+
+### 2. **Early Warning System**
+   - **Anomaly Detection:** Identify unexpected spikes in heavy metal levels to prevent potential regulatory violations or equipment failures.
+   - **Compliance Monitoring:** Predict when heavy metal concentrations approach regulatory limits and take preventive action.
+
+### 3. **Optimizing Chemical Usage**
+   - **Predictive Dosing Optimization:** Forecast the necessary chemical dosage to maintain target heavy metal levels, reducing chemical costs.
+   - **Focus on Key Influencers:** Utilize sensitivity analysis results to prioritize influential parameters (e.g., electrical conductivity) for better control.
+
+### 4. **Scenario Analysis**
+   - **Simulate Operating Conditions:** Evaluate the impact of different treatment strategies and operating conditions on heavy metal removal efficiency.
+   - **Assess Upgrades:** Predict the effects of potential system upgrades or changes in the treatment process.
+
+### 5. **Deployment Considerations**
+   - **Model Retraining:** Periodically update the model with new data to maintain accuracy.
+   - **Cloud vs. Edge Deployment:** Choose deployment based on latency requirements—cloud for centralized processing or edge for faster local predictions.
+
+These recommendations help guide the practical use of the PSO-LSTM model for optimizing wastewater treatment and ensuring compliance with regulatory standards.
+
 ## License
 This project is licensed under the MIT License - see the LICENSE file for details.
 
diff --git a/var.py b/var.py
new file mode 100644
index 0000000..f7ccf5d
--- /dev/null
+++ b/var.py
@@ -0,0 +1,83 @@
+# -*- coding: utf-8 -*-
+"""
+VAR Model for Time Series Forecasting
+
+This script performs Vector Autoregression (VAR) for time series forecasting
+of multiple variables related to industrial wastewater treatment. 
+It includes data loading, model fitting, and predictions.
+"""
+
+# Import required packages
+import pandas as pd
+import matplotlib.pyplot as plt
+from statsmodels.tsa.vector_ar.var_model import VAR
+
+# Load the dataset
+dataset = pd.read_excel("202106有機自動加藥數據統計(新).xlsx", header=0, index_col=0)
+# Select relevant columns
+dataset = dataset[['銅濃(mg/L)', '銅在線濃度(mg/L)', 'OOO', 'pH', 'ph槽ORP', 'Chemical_A_ORP', 'Chemical_A', 'Chemical_B']]
+
+# Display the first few rows of the dataset
+print(dataset.head())
+
+# Check data types for each column
+print("\nData Types:\n", dataset.dtypes)
+
+# Split the data into training and validation sets (80% train, 20% validation)
+train_size = int(0.8 * len(dataset))
+train = dataset[:train_size]
+valid = dataset[train_size:]
+
+# Fit the VAR model on the training set
+model = VAR(endog=train)
+model_fit = model.fit()
+
+# Summary of the model
+print("\nModel Summary:\n")
+print(model_fit.summary())
+
+# Make predictions on the validation set
+predictions = model_fit.forecast(y=model_fit.y, steps=len(valid))
+
+# Convert predictions to DataFrame for better readability
+predicted_df = pd.DataFrame(predictions, index=valid.index, columns=valid.columns)
+
+# Plotting the actual vs predicted values for the validation set
+plt.figure(figsize=(12, 6))
+for column in valid.columns:
+    plt.plot(valid.index, valid[column], label=f'Actual {column}')
+    plt.plot(predicted_df.index, predicted_df[column], linestyle='--', label=f'Predicted {column}')
+plt.title('Actual vs Predicted Values for Validation Set')
+plt.xlabel('Time')
+plt.ylabel('Value')
+plt.legend(loc='upper left', bbox_to_anchor=(1, 1))
+plt.grid()
+plt.tight_layout()
+plt.show()
+
+# Make final predictions on the entire dataset
+final_model = VAR(endog=dataset)
+final_model_fit = final_model.fit()
+future_steps = 10  # Number of future steps to predict
+final_predictions = final_model_fit.forecast(y=final_model_fit.y, steps=future_steps)
+
+# Convert the final predictions to a DataFrame
+future_index = pd.date_range(start=dataset.index[-1], periods=future_steps + 1, freq='D')[1:]
+final_predicted_df = pd.DataFrame(final_predictions, index=future_index, columns=dataset.columns)
+
+# Display the final forecasted values
+print("\nFinal Forecasted Values:\n")
+print(final_predicted_df)
+
+# Plotting the future predictions
+plt.figure(figsize=(12, 6))
+for column in dataset.columns:
+    plt.plot(dataset.index, dataset[column], label=f'Historical {column}')
+    plt.plot(final_predicted_df.index, final_predicted_df[column], linestyle='--', label=f'Forecasted {column}')
+plt.title('Future Forecasted Values')
+plt.xlabel('Time')
+plt.ylabel('Value')
+plt.legend(loc='upper left', bbox_to_anchor=(1, 1))
+plt.grid()
+plt.tight_layout()
+plt.show()