Homework03

Introduction

The following RMD contains CUNY SPS DATA 605 Fall 2025 Homework 03 assignment. This report contains five problems.

Load Libraries

library(tidyverse)
library(knitr)
library(nortest)
library(dplyr)
library(glue)
library(gt)
library(lmtest)
library(car) 
library(broom)

Problem 1: Transportation Safety

# Peak at the cars dataset
head(cars)

##   speed dist
## 1     4    2
## 2     4   10
## 3     7    4
## 4     7   22
## 5     8   16
## 6     9   10

summary(cars)

##      speed           dist       
##  Min.   : 4.0   Min.   :  2.00  
##  1st Qu.:12.0   1st Qu.: 26.00  
##  Median :15.0   Median : 36.00  
##  Mean   :15.4   Mean   : 42.98  
##  3rd Qu.:19.0   3rd Qu.: 56.00  
##  Max.   :25.0   Max.   :120.00

The cars dataset gives the speed of cars and the distances taken to stop and has data recorded in the 1920s. The data frame shows two fields, speed and distance (speed and dist). Speed ranges from 4mph to 25mph and distance ranges from 2ft to 120ft

Question 1: Data Visualization:

# Create the plot with regression line
ggplot(cars, aes(x = speed, y = dist, alpha = 0.7)) +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE, color = "salmon") +
  labs(title = "Car Stopping Distance vs. Speed",
       x = "Speed (mph)",
       y = "Stopping Distance (ft)") +
  theme_minimal()

The trend here looks curved, meaning the relationship between the two may not be perfectly linear. There is also an increasing trend seen.

Question 2: Build a Linear Model:

# Build a basic linear regression model
car_model <- lm(dist ~ speed, data = cars, na.action = na.exclude)
model_summary <- summary(car_model)
model_summary

## 
## Call:
## lm(formula = dist ~ speed, data = cars, na.action = na.exclude)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -29.069  -9.525  -2.272   9.215  43.201 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -17.5791     6.7584  -2.601   0.0123 *  
## speed         3.9324     0.4155   9.464 1.49e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.38 on 48 degrees of freedom
## Multiple R-squared:  0.6511, Adjusted R-squared:  0.6438 
## F-statistic: 89.57 on 1 and 48 DF,  p-value: 1.49e-12

The model’s p-value (far less than 0.001) shows that speed is quite the statistically significant predictor of stopping distance.

Question 3: Model Quality Evaluation:

# Get R-squared 
r_squared <- model_summary$r.squared
cat("R-squared value:", r_squared, "\n")

## R-squared value: 0.6510794

This R-squared value means that 65.11% of the variance in stopping distance is explained by speed. This R-squared is moderate, not incredible. The remaining 34.9% of variance is likely due to factors not included in this data set like driver behavior, braking system differences, road surface conditions, and measurement variation-factors.

Question 4: Residual Analysis:

# Plot residuals vs fitted
plot(car_model$fitted.values, car_model$residuals,
     xlab = "Fitted Values", ylab = "Residuals",
     main = "Residuals vs Fitted Values")
abline(h = 0, col = "salmon", lty = 2)

The funnel shape suggests heteroscedasticity. The curve suggests non-linearity.

# Create a Q-Q plot for residuals
qqnorm(car_model$residuals)
qqline(car_model$residuals, col = "salmon")

The tails deviate on both sides, most severely on the right of the graph. This suggests some right skewedness.

# Anderson-Darling test
ad_test <- ad.test(car_model$residuals)
print(ad_test)

## 
##  Anderson-Darling normality test
## 
## data:  car_model$residuals
## A = 0.79406, p-value = 0.0369

Because the p-value is below 0.05, the null hypothesis that the residuals are normally distributed is rejected. This test agrees with the Q-Q plot’s findings of skewedness.

# Histogram of residuals
hist(car_model$residuals, breaks = 15, 
     main = "Histogram of Residuals",
     xlab = "Residuals", 
     col = "pink",
     border = "gray")

The histogram shows a visualization of the right-skewedness, and has a gap in the middle where the Q-Q plot shows a jump in values as well.

Question 5: Conclusion:

Based on the findings above, the linear model does capture the general upward trend but fails to account for curvature and increasing variance at higher speeds. The model explains ~65% of variance (R-squared = 0.651), but residual diagnostics suggest non-linearity and heteroscedasticity (AD p = 0.0369; funnel-shaped residuals). This suggests that the physical braking process is not perfectly linear-consistent with real world stopping physics, where kinetic energy increases with the square of speed. A polynomial/quadratic model is likely more appropriate.

Here is a polynomial model for comparison:

# Fit polynomial model
car_model2 <- lm(dist ~ speed + I(speed^2), data = cars)
summary(car_model2)

## 
## Call:
## lm(formula = dist ~ speed + I(speed^2), data = cars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -28.720  -9.184  -3.188   4.628  45.152 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  2.47014   14.81716   0.167    0.868
## speed        0.91329    2.03422   0.449    0.656
## I(speed^2)   0.09996    0.06597   1.515    0.136
## 
## Residual standard error: 15.18 on 47 degrees of freedom
## Multiple R-squared:  0.6673, Adjusted R-squared:  0.6532 
## F-statistic: 47.14 on 2 and 47 DF,  p-value: 5.852e-12

# Get smooth prediction line and CI
new_speed <- data.frame(speed = seq(min(cars$speed), max(cars$speed), length.out = 200))
preds <- predict(car_model2, newdata = new_speed, se.fit = TRUE)
new_speed$fit <- preds$fit
new_speed$upper <- preds$fit + 1.96 * preds$se.fit
new_speed$lower <- preds$fit - 1.96 * preds$se.fit

# Vizualise with linear fit
ggplot(cars, aes(x = speed, y = dist)) +
  geom_point(alpha = 0.7) +
  geom_ribbon(data = new_speed, aes(x = speed, ymin = lower, ymax = upper),
              inherit.aes = FALSE, alpha = 0.15) +
  geom_line(data = new_speed, 
            aes(x = speed, y = fit), 
            inherit.aes = FALSE,
            linewidth = 1, linetype = "solid",
            color = "aquamarine3"
            ) +
  geom_smooth(method = "lm", formula = y ~ x, se = FALSE, color = "salmon",
              linetype = "dashed", size = 0.9) +
  labs(title = "Stopping Distance vs Speed — Quadratic fit (blue) vs Linear fit (dashed salmon)",
       subtitle = "Blue line = quadratic model with 95% CI; dashed salmon = simple linear model",
       x = "Speed (mph)",
       y = "Stopping distance (ft)") +
  theme_minimal()

R-squared increased slightly to 0.667; this suggests curvature is present. However, the residual SD changed only marginally. The next step would be to check influence points and compare AIC/BIC.

Problem 2: Health Policy Analyst

# Download provided data set
who_raw <- read_csv("https://github.com/evanskaylie/DATA605/raw/refs/heads/main/who.csv")

# Save the data 
who_data <- who_raw

# Take a look
head(who_data)

## # A tibble: 6 × 10
##   Country   LifeExp InfantSurvival Under5Survival TBFree  PropMD  PropRN PersExp
##   <chr>       <dbl>          <dbl>          <dbl>  <dbl>   <dbl>   <dbl>   <dbl>
## 1 Afghanis…      42          0.835          0.743  0.998 2.29e-4 5.72e-4      20
## 2 Albania        71          0.985          0.983  1.00  1.14e-3 4.61e-3     169
## 3 Algeria        71          0.967          0.962  0.999 1.06e-3 2.09e-3     108
## 4 Andorra        82          0.997          0.996  1.00  3.30e-3 3.5 e-3    2589
## 5 Angola         41          0.846          0.74   0.997 7.04e-5 1.15e-3      36
## 6 Antigua …      73          0.99           0.989  1.00  1.43e-4 2.77e-3     503
## # ℹ 2 more variables: GovtExp <dbl>, TotExp <dbl>

str(who_data)

## spc_tbl_ [190 × 10] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ Country       : chr [1:190] "Afghanistan" "Albania" "Algeria" "Andorra" ...
##  $ LifeExp       : num [1:190] 42 71 71 82 41 73 75 69 82 80 ...
##  $ InfantSurvival: num [1:190] 0.835 0.985 0.967 0.997 0.846 0.99 0.986 0.979 0.995 0.996 ...
##  $ Under5Survival: num [1:190] 0.743 0.983 0.962 0.996 0.74 0.989 0.983 0.976 0.994 0.996 ...
##  $ TBFree        : num [1:190] 0.998 1 0.999 1 0.997 ...
##  $ PropMD        : num [1:190] 2.29e-04 1.14e-03 1.06e-03 3.30e-03 7.04e-05 ...
##  $ PropRN        : num [1:190] 0.000572 0.004614 0.002091 0.0035 0.001146 ...
##  $ PersExp       : num [1:190] 20 169 108 2589 36 ...
##  $ GovtExp       : num [1:190] 92 3128 5184 169725 1620 ...
##  $ TotExp        : num [1:190] 112 3297 5292 172314 1656 ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   Country = col_character(),
##   ..   LifeExp = col_double(),
##   ..   InfantSurvival = col_double(),
##   ..   Under5Survival = col_double(),
##   ..   TBFree = col_double(),
##   ..   PropMD = col_double(),
##   ..   PropRN = col_double(),
##   ..   PersExp = col_double(),
##   ..   GovtExp = col_double(),
##   ..   TotExp = col_double()
##   .. )
##  - attr(*, "problems")=<externalptr>

TotExp is total health expenditure in USD; LifeExp is years. Note wide range and potential outliers (very large TotExp)

Question 1: Initial Assessment of Healthcare Expenditures and Life Expectancy

# Create plot for expenditures vs life expectancy
ggplot(who_data, mapping =  aes(x = TotExp, y = LifeExp)) +
  geom_point(alpha = 0.7) +
  theme_minimal() +
  labs(title = "Scatterplot of Life Expectancy vs Life ExpectancyTotal Healthcare Expenditure",
       x = "Total Healthcare Expenditure (TotExp)",
       y = "Life Expectancy (LifeExp)") +
  geom_smooth(method = "lm", se = FALSE, col = "salmon") +
  theme_minimal()

Visually, this does not look like a linear relationship at all. The variance increases as life expectancy and expenditure increase, as well as a dare-I-say, quadratic relationship.

# Run simple linear regression
who_model1 <- lm(LifeExp ~ TotExp, data = who_data)
who_model1_summary <- summary(who_model1)

# Tidy table for future use
who_model1_tidy <- tidy(who_model1)
who_model1_glance <- glance(who_model1)

# Display model summary
who_model1_summary

## 
## Call:
## lm(formula = LifeExp ~ TotExp, data = who_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -24.764  -4.778   3.154   7.116  13.292 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 6.475e+01  7.535e-01  85.933  < 2e-16 ***
## TotExp      6.297e-05  7.795e-06   8.079 7.71e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.371 on 188 degrees of freedom
## Multiple R-squared:  0.2577, Adjusted R-squared:  0.2537 
## F-statistic: 65.26 on 1 and 188 DF,  p-value: 7.714e-14

# Extract key statistics from data frame in a table
who_stats_table1 <- data.frame(
  Statistic = c("F-statistic", "F-test p-value", "R-squared", "Residual Std. Error"),
  Value = c(
    who_model1_summary$fstatistic[1],
    pf(who_model1_summary$fstatistic[1],
       who_model1_summary$fstatistic[2],
       who_model1_summary$fstatistic[3],
       lower.tail = FALSE),
    who_model1_summary$r.squared,
    who_model1_summary$sigma
  )
)

# Print as a formatted table
kable(who_stats_table1, caption = "Key Statistics for Linear Regression Model")

Key Statistics for Linear Regression Model

Statistic	Value
F-statistic	65.2641982
F-test p-value	0.0000000
R-squared	0.2576922
Residual Std. Error	9.3710333

# Check regression assumptions
par(mfrow = c(2, 2))
plot(who_model1)

# Anderson-Darling test
who_ad_test <- ad.test(who_model1$residuals)
print(who_ad_test)

## 
##  Anderson-Darling normality test
## 
## data:  who_model1$residuals
## A = 7.2968, p-value < 2.2e-16

The AD test p-value (r round(who_ad_test$p.value, 4)) suggests residuals deviate from normality.

# Breusch-Pagan test for heteroscedasticity
bptest(who_model1)

## 
##  studentized Breusch-Pagan test
## 
## data:  who_model1
## BP = 2.6239, df = 1, p-value = 0.1053

The Breusch–Pagan test returned p = 0.1053 so we do not reject constant variance at α = 0.05. However visual diagnostics and the AD test show non-normal residuals and possible spread at high TotExp; combined evidence suggests we should explore transformations and influence points.

# Influence: Cooks distance
cd1 <- cooks.distance(who_model1)
plot(cd1, type = "h", main = "Cook's distance (model1)")
abline(h = 4/(nrow(who_data)-length(coef(who_model1))), col = "salmon")

which(cd1 > 4/(nrow(who_data)-length(coef(who_model1))))

##  76  99 111 125 152 
##  76  99 111 125 152

There seem to be a series of outliers affecting the data here.

Question 2: Transforming Variables for a Better Fit

# Transform variables
who_data$LifeExp_trans <- who_data$LifeExp^4.6
who_data$TotExp_trans <- who_data$TotExp^0.06

# Create scatterplot with transformed variables
ggplot(who_data, aes(x = TotExp_trans, y = LifeExp_trans)) +
  geom_point(alpha = 0.7) +
  geom_smooth(method = "lm", se = TRUE, color = "salmon") +
  labs(title = "Transformed Total Healthcare Expenditure vs Life Expectancy",
       x = "Transformed Life Expectancy (LifeExp^4.6)",
       y = "Transformed Total Healthcare Expenditure (TotExp^0.06)") +
  theme_minimal()

In this step, power transforms are applied to linearize curvature. A warning here: small changes in transformed space map to extreme changes in original TotExp. Interpret back-transformed forecasts only within the observed TotExp range. For extrapolations use prediction intervals via bootstrapping.

# Run linear regression with transformed variables
who_model2 <- lm(LifeExp_trans ~ TotExp_trans, data = who_data)
who_model2_summary <- summary(who_model2)

# Display model summary
who_model2_summary

## 
## Call:
## lm(formula = LifeExp_trans ~ TotExp_trans, data = who_data)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -308616089  -53978977   13697187   59139231  211951764 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -736527910   46817945  -15.73   <2e-16 ***
## TotExp_trans  620060216   27518940   22.53   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 90490000 on 188 degrees of freedom
## Multiple R-squared:  0.7298, Adjusted R-squared:  0.7283 
## F-statistic: 507.7 on 1 and 188 DF,  p-value: < 2.2e-16

# Extract key statistics from data frame to a table
who_stats_table2 <- data.frame(
  Statistic = c("F-statistic", "F-test p-value", "R-squared", "Residual Std. Error"),
  Value = c(
    who_model2_summary$fstatistic[1],
    pf(who_model2_summary$fstatistic[1],
       who_model2_summary$fstatistic[2],
       who_model2_summary$fstatistic[3],
       lower.tail = FALSE),
    who_model2_summary$r.squared,
    who_model2_summary$sigma
  )
)

# Print as a formatted table
kable(who_stats_table2, caption = "Key Statistics for Polynomial Regression Model")

Key Statistics for Polynomial Regression Model

Statistic	Value
F-statistic	5.076967e+02
F-test p-value	0.000000e+00
R-squared	7.297673e-01
Residual Std. Error	9.049239e+07

# Check regression assumptions for transformed model
par(mfrow = c(2, 2))
plot(who_model2)

# Anderson-Darling test
who_ad_test2 <- ad.test(who_model2$residuals)
print(who_ad_test2)

## 
##  Anderson-Darling normality test
## 
## data:  who_model2$residuals
## A = 1.4003, p-value = 0.001236

# Question 1 linear results
kable(who_stats_table1, caption = "Key Statistics for Linear Regression Model")

Key Statistics for Linear Regression Model

Statistic	Value
F-statistic	65.2641982
F-test p-value	0.0000000
R-squared	0.2576922
Residual Std. Error	9.3710333

# Question 2 transfomed results
kable(who_stats_table2, caption = "Key Statistics for Polynomial Regression Model")

Key Statistics for Polynomial Regression Model

Statistic	Value
F-statistic	5.076967e+02
F-test p-value	0.000000e+00
R-squared	7.297673e-01
Residual Std. Error	9.049239e+07

Model 1 shows strong non-normality (AD p < 0.001) and heteroscedasticity, with residuals spreading at higher TotExp. Model 2’s transformation reduces these issues: residuals are closer to normal (AD p ≈ 0.001) and variance is more uniform, though a few influential points remain. Overall, the transformation improves model assumptions and fit.

R-squared increased in the transformed model because the power transformations better captured the non-linear relationship between TotExp and LifeExp. While the model is highly significant statistically, large effect sizes may not always translate to practically meaningful life expectancy gains for modest expenditure changes.

Question 3: Forecasting Life Expectancy Based on Transformed Expenditures

# Extract model coefficients
who_intercept2 <- who_model2$coefficients[1]
who_slope2 <- who_model2$coefficients[2]

# Forecast life expectancy for TotExp^0.06 = 1.5
who_forecast1_trans <- who_intercept2 + who_slope2 * 1.5
who_forecast1 <- who_forecast1_trans^(1/4.6)

# Forecast life expectancy for TotExp^0.06 = 2.5
who_forecast2_trans <- who_intercept2 + who_slope2 * 2.5
who_forecast2 <- who_forecast2_trans^(1/4.6)

# Results
cat("Forecasted life expectancy when TotExp^0.06 = 1.5:", who_forecast1, "years\n")

## Forecasted life expectancy when TotExp^0.06 = 1.5: 63.31153 years

cat("Forecasted life expectancy when TotExp^0.06 = 2.5:", who_forecast2, "years\n")

## Forecasted life expectancy when TotExp^0.06 = 2.5: 86.50645 years

# Calculate original TotExp values for these forecasts
who_totexp1 <- 1.5^(1/0.06)
who_totexp2 <- 2.5^(1/0.06)
cat("This corresponds to TotExp =", who_totexp1, "USD for forecast 1\n")

## This corresponds to TotExp = 860.705 USD for forecast 1

cat("This corresponds to TotExp =", who_totexp2, "USD for forecast 2\n")

## This corresponds to TotExp = 4288777 USD for forecast 2

Back-transforming TotExp^0.06 can produce extreme USD values. Interpret forecasts only within realistic expenditure ranges

The forecasts from the transformed model show that healthcare spending is strongly linked to life expectancy, but the relationship is non-linear. This means that increases in spending do not always lead to proportional gains in life expectancy. The model helps guide policy by showing how the impact of spending can vary depending on a country’s current level of healthcare investment.

Question 4: Interaction Effects in Multiple Regression

# Create interaction term
who_data$PropMD_TotExp <- who_data$PropMD * who_data$TotExp

# Build multiple regression model with interaction
who_model3 <- lm(LifeExp ~ PropMD + TotExp + PropMD_TotExp, data = who_data)
who_model3_summary <- summary(who_model3)

# Display model summary
who_model3_summary

## 
## Call:
## lm(formula = LifeExp ~ PropMD + TotExp + PropMD_TotExp, data = who_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -27.320  -4.132   2.098   6.540  13.074 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    6.277e+01  7.956e-01  78.899  < 2e-16 ***
## PropMD         1.497e+03  2.788e+02   5.371 2.32e-07 ***
## TotExp         7.233e-05  8.982e-06   8.053 9.39e-14 ***
## PropMD_TotExp -6.026e-03  1.472e-03  -4.093 6.35e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.765 on 186 degrees of freedom
## Multiple R-squared:  0.3574, Adjusted R-squared:  0.3471 
## F-statistic: 34.49 on 3 and 186 DF,  p-value: < 2.2e-16

# Extract key statistics from data frame to a table
who_stats_table3 <- data.frame(
  Statistic = c("F-statistic", "F-test p-value", "R-squared", "Residual Std. Error"),
  Value = c(
    who_model3_summary$fstatistic[1],
    pf(who_model3_summary$fstatistic[1],
       who_model3_summary$fstatistic[2],
       who_model3_summary$fstatistic[3],
       lower.tail = FALSE),
    who_model3_summary$r.squared,
    who_model3_summary$sigma
  )
)

# Print as a formatted table
kable(who_stats_table3, caption = "Key Statistics for Multiple Regression Model")

Key Statistics for Multiple Regression Model

Statistic	Value
F-statistic	34.4883268
F-test p-value	0.0000000
R-squared	0.3574352
Residual Std. Error	8.7654934

The negative interaction means that the effect of healthcare spending on life expectancy is smaller in countries with more doctors. This shows that the impact of resources isn’t straightforward and depends on existing workforce levels. The recommendation is to take a balanced approach, investing in both healthcare spending and staffing to most effectively improve population health.

Question 5: Forecasting Life Expectancy with Interaction Terms

# Extract model coefficients
# Extract coefficients from model
who_beta0 <- who_model3$coefficients["(Intercept)"]
who_beta1 <- who_model3$coefficients["PropMD"]
who_beta2 <- who_model3$coefficients["TotExp"]
who_beta3 <- who_model3$coefficients["PropMD_TotExp"]

# Forecasting values
who_propmd_val <- 0.03
who_totexp_val <- 14

# Forecast life expectancy
who_predicted_lifeexp <- who_beta0 + 
                         who_beta1 * who_propmd_val + 
                         who_beta2 * who_totexp_val + 
                         who_beta3 * (who_propmd_val * who_totexp_val)

# Print results
cat("Predicted life expectancy for a country with:\n")

## Predicted life expectancy for a country with:

cat("PropMD =", who_propmd_val, "\n")

## PropMD = 0.03

cat("TotExp =", who_totexp_val, "\n")

## TotExp = 14

cat("Predicted Life Expectancy:", round(who_predicted_lifeexp, 1), "years\n")

## Predicted Life Expectancy: 107.7 years

# Compare to average in dataset
who_avg_lifeexp <- mean(who_data$LifeExp, na.rm = TRUE)
cat("Average life expectancy in dataset:", round(who_avg_lifeexp, 1), "years\n")

## Average life expectancy in dataset: 67.4 years

Life expectancy increases with both more doctors and higher spending. The negative interaction means the benefit of extra spending decreases when doctor density is high (and vice versa). Predictions should be interpreted cautiously if the values are outside the range observed in the dataset. The predicted life expectency is 107.7 years. I would say this is a bit unrealistic and may be caused by extrapolation.

Problem 3: Retail Company Analyst

Question 1: Inventory Cost

Get derivative and set it to zero

The given cost function is:

\[C(Q)=\frac{D}{Q}S + \frac{Q}{2}H\]

Differentiate with respect to \(Q\) and set to zero:

\[C'(Q) = -\frac{D S}{Q^2} + \frac{H}{2} = 0\]

Solve for the optimal lot size \(Q^*\):

\[C'(Q)= -\frac{DS}{Q^2} + \frac{H}{2}=0\] \[\frac{H}{2} = \frac{DS}{Q^2}\]

\[\frac{H}{2}Q^2 = DS\]

\[HQ^2 = 2DS\]

\[Q^2 = \frac{2DS}{H}\]

\[Q^* = \sqrt{\frac{2DS}{H}}\]

Provided parameters

• Annual demand: D = 110 units
• Ordering cost: S = $8.25 per order
• Holding cost: H = $3.75 per unit/year

Final answers (plug in values)

• Optimal order quantity Q* = 22.000 units
• Orders per year N = 5.000
• Total annual inventory cost at Q* = $82.50

In summary, the optimal lot size is 22 units and the number of orders the company should place per year to minimize total inventory cost is 5.

Question 2: Revenue Maximization

# Given revenue function
R <- function(t) { -3150 * t^(-4) - 220 * t + 6530 }

\[12600 t^{-5} - 220 = 0\] \[12600 / t^5 = 220\] \[t^5 = 12600 / 220\] \[t^* = (630/11)^{1/5} pprox 2.2469\]

# Save critical point
t_star <- (630/11)^(1/5)

# Do the evaluation
t_star_approx <- unname(t_star)
R_max <- R(t_star_approx)

# Print results
cat(sprintf("Critical time: t* ≈ %.6f days\n", t_star_approx))

Critical time: t* ≈ 2.246930 days

cat(sprintf("Maximum revenue: R(t*) ≈ $%.2f\n", R_max))

Maximum revenue: R(t*) ≈ $5912.09

# Now verify this is the maximum by checking the second derivitive here
R2 <- function(t) { -63000 * t^(-6) }
cat(sprintf("Second derivative at t*: R''(t*) ≈ %.4e (negative => maximum)\n", R2(t_star_approx)))

## Second derivative at t*: R''(t*) ≈ -4.8956e+02 (negative => maximum)

Because this R’’(t*) is less than 0, the critical point above is a maximum.

# Now plug the 2.25 value for t* in to R(t)
R_at_t_star <- R(t_star_approx)
cat(sprintf("Revenue at t = 2.25: R(2.25) = $%.2f\n", R_at_t_star))

## Revenue at t = 2.25: R(2.25) = $5912.09

Here, this tells us the maximum revenue the company can expect from this campaign is about $5912.09.

Question 3: Demand Area Under Curve

Total Revenue from Demand Curve

Given demand function:

\[P(x) = 2x - 9.3\]

To find total revenue between (x_1 = 2) and (x_2 = 5), compute the definite integral:

\[\text{Revenue} = \int_{2}^{5} (2x - 9.3) \, dx\]

Evaluating the integral:

\[[x^2 - 9.3x]_{2}^{5}\]

Total revenue over x = 2 to x = 5:

$-6.90 dollars

This revenue is negative, because this function has integration over negative prices. Negative revenue is not ideal, so the suggestion would be to not sell the product between the two given quantity levels, 2 and 5. The product is more expensive to supply than the demand is willing to pay, so there is revenue loss.

Question 4: Profit Optimization

# Given profit function
pi_fun <- function(x) { x * log(9 * x) - (x^6) / 6 }
# Pi'(x) to set to zero
pi_prime <- function(x) { log(9 * x) + 1 - x^5 } 

# Next, set pi_prime to zero to find the critical point 

## This line below actually does not run, because the function does not cross zero

# crit_point <- uniroot(pi_prime, lower = 0.01, upper = 2)$root

# Instead, visualize the function 
curve(log(9*x) + 1 - x^5, from = 0.01, to = 5, ylab = "pi'(x)")
abline(h = 0, col = "salmon")

The profit function is decreasing on the feasible domain and no positive x yields a maximum.

Question 5: Spending Behavior

# Given probability density function 
f_q5 <- function(x) { 1 / (6 * x) }

# Given interval parameters
a <- 1
b <- exp(6)

# If integral over this range equals one, it is a valid probability density function
total_prob <- integrate(f_q5, lower = a, upper = b)$value
# Conditional output
if (abs(total_prob - 1) < 1e-8) {
  cat("The function f(x) = 1/(6x) is a valid probability density function.\n")
} else {
  cat("The function f(x) = 1/(6x) is NOT a valid probability density function.\n")
}

## The function f(x) = 1/(6x) is a valid probability density function.

# Print the probability over the interval
cat(sprintf("Total probability over [1, e^6] = %.2f\n", total_prob))

## Total probability over [1, e^6] = 1.00

# Now calculate the probability that customer spends between 1 and e^6
prob_interval <- integrate(f_q5, lower = a, upper = b)$value
cat(sprintf("Probability a customer spends between 1 and e^6 is %.2f\n", prob_interval))

## Probability a customer spends between 1 and e^6 is 1.00

This spending probability makes sense as the spending range is the full range of the probability density function.

Question 6: Market Share Estimation

# Given rate of market penetration
dN_dt <- function(t) { 500 / (t^4 + 10) }

# Initial condition given parameters
t0 <- 1
N0 <- 6530

# Function to compute cumulative market share N(t)
N_t <- function(t) {
  # Integrate dN/dt from t0 to t
  delta_N <- integrate(dN_dt, lower = t0, upper = t)$value
  # Add initial market share
  N0 + delta_N
}

# Calculate market share after 10 days
t_final <- 10
market_share_10 <- N_t(t_final)

# Print results
cat(sprintf("Cumulative market share at t = %d days: %.2f units\n", t_final, market_share_10))

## Cumulative market share at t = 10 days: 6579.54 units

Problem 4: Business Optimization

Question 1: Revenue and Cost

Approximating the Revenue Function

# Given revenue funcition
cat("Revenue function: R(x) = e^x\n")

## Revenue function: R(x) = e^x

# Maclaurin series coefficients for e^x
R0 <- exp(0) # R(0)
R1 <- exp(0) # R'(0)
R2 <- exp(0)/2 # R''(0)/2 factorial for second-degree
cat(sprintf("R(0) = %.0f, R'(0) = %.0f, R''(0)/2 = %.1f\n", R0, R1, R2))

## R(0) = 1, R'(0) = 1, R''(0)/2 = 0.5

# Maclaurin series approximation
R_approx <- function(x) { R0 + R1*x + R2*x^2 }
cat("Maclaurin series approximation: R(x) ≈ 1 + x + 0.5*x^2\n\n")

## Maclaurin series approximation: R(x) ≈ 1 + x + 0.5*x^2

The Maclaurin Series gives a simple approximation of revenue for small sales volumes. Using R(x) ≈ 1 + x + 0.5*x^2, I estimate revenue without computing the exponential function. This is helpful for rapid “what-if” scenarios in sales planning.

Approximating the Cost Function

# Given cost function
cat("Cost function: C(x) = ln(1 + x)\n")

## Cost function: C(x) = ln(1 + x)

# Maclaurin series coefficients for ln(1+x)
C0 <- log(1 + 0) # C(0)
C1 <- 1/(1 + 0) # C'(0)
C2 <- -1/(2*(1 + 0)^2) # C''(0)/2
cat(sprintf("C(0) = %.0f, C'(0) = %.0f, C''(0)/2 = %.2f\n", C0, C1, C2))

## C(0) = 0, C'(0) = 1, C''(0)/2 = -0.50

# Maclaurin series approximation
C_approx <- function(x) { C0 + C1*x + C2*x^2 }
cat("Maclaurin series approximation: C(x) ≈ 0 + x - 0.5x^2 = x - 0.5x^2\n\n")

## Maclaurin series approximation: C(x) ≈ 0 + x - 0.5x^2 = x - 0.5x^2

The Maclaurin Series approximates production cost for small quantities. Using C(x) ≈ x - 0.5*x^2 simplifies calculations and helps managers quickly estimate cost changes as production increases. It highlights that initial cost increases almost linearly with small production.

Linear vs Nonlinear Optimization

# Original profit function
profit_original <- function(x) { exp(x) - log(1 + x) }

# Approximated profit function
profit_approx <- function(x) { R_approx(x) - C_approx(x) }
cat("Approximated profit function: PI_approx(x) = 1 + x^2\n\n")

## Approximated profit function: PI_approx(x) = 1 + x^2

The second-degree Taylor approximation suggests PI_approx(x) = 1 + x-squared, which is increasing in x for x > 0, so the approximation does not produce a sensible finite interior maximizer. This highlights that polynomial approximations can change optimization structure; always compare with the original function numerically.

# Define approximated revenue and cost functions
R_approx <- function(x) { 1 + x + 0.5*x^2 }
C_approx <- function(x) { x - 0.5*x^2 }

# Approximate profit
profit_approx <- function(x) { R_approx(x) - C_approx(x) }

# Analytical derivative of approximated profit
profit_approx_prime <- function(x) { 2*x }

# Find maximum of original profit numerically
max_profit <- optimize(profit_original, interval=c(0,2), maximum=TRUE)

# Solve derivative = 0
x_opt_approx <- 0
profit_opt_approx <- profit_approx(x_opt_approx)

# Results
cat("Approximate optimal production:", x_opt_approx, "\n")

## Approximate optimal production: 0

cat("Profit at approximate optimum:", profit_opt_approx, "\n")

## Profit at approximate optimum: 1

x_vals <- seq(0, 2, by=0.01)

df_p4 <- data.frame(
                  x = x_vals,
                  Original = sapply(x_vals, profit_original),
                  Approximation = sapply(x_vals, profit_approx)
                  )

# Visualize
ggplot(df_p4, aes(x=x)) +
  geom_line(aes(y=Original, color="Original")) +
  geom_line(aes(y=Approximation, color="Approximation")) +
  geom_vline(xintercept = max_profit$maximum, linetype="dashed", color="aquamarine3") +
  geom_vline(xintercept = x_opt_approx, linetype="dotted", color="salmon") +
  labs(title="Profit: Original vs Approximated Function",
      x="Units Sold",
      y="Profit",
      color="Function") +
  theme_minimal()

The plot compares original vs approximated profit. The vertical lines mark optimal production for each method. Managers can quickly see where small approximations give similar results to the full calculation and understand the safe range to rely on simpler estimates.

Question 2: Financial Modeling

Maclaurin Series Expansion

# Original risk function
risk <- function(x) sqrt(x)

# Note: sqrt(x) is not differentiable at x = 0
# We'll use a small positive x0 for linear approximation
x0 <- 0.01

# Linear (first-order Taylor) approximation around x0
risk_approx <- function(x) {
  risk(x0) + (1/(2*sqrt(x0))) * (x - x0)
}

risk_approx(0.01)  # Check approximation at x0

## [1] 0.1

The square root function cannot be expanded exactly at 0 because derivatives are undefined. For very small investments, we can use a linear approximation around a small x0, which is easy to calculate. Because sqrt(x) is not differentiable at zero, we approximate around a small x0. This gives a simple linear formula to estimate risk for very small investments, helping analysts make fast decisions under low-risk scenarios. Quick note: sqrt(x) has derivative infinite at 0 so Maclaurin at 0 is not defined; linearizing around a small x0 gives a local approximation only.

Practical Application

# Investment amounts to check
x_vals <- c(0.01, 0.05, 0.1, 0.5, 1)

# Compare actual vs approximated risk
risk_df <- data.frame(
Investment = x_vals,
Actual_Risk = risk(x_vals),
Approximated_Risk = risk_approx(x_vals)
)

risk_df

##   Investment Actual_Risk Approximated_Risk
## 1       0.01   0.1000000              0.10
## 2       0.05   0.2236068              0.30
## 3       0.10   0.3162278              0.55
## 4       0.50   0.7071068              2.55
## 5       1.00   1.0000000              5.05

The table shows that for very small investments (x < 0.1), the linear approximation is very close to the real risk. For moderate investments, the approximation diverges. Managers can use this to understand the limits of quick estimates.

Optimization Scenario

# Set maximum risk allowed
risk_tolerance <- 0.2

# Solve for optimal investment using linear approximation
x_opt <- x0 + (risk_tolerance - risk(x0)) / (1/(2*sqrt(x0)))
x_opt

## [1] 0.03

This calculation gives the largest safe investment under a given risk tolerance. Using the linear approximation is fast and sufficient for small investments, but for larger amounts, the exact function should be used to avoid underestimating risk.

Visualization

# Generate sequence for plotting
x_plot <- seq(0.001, 1, by=0.001)

# Plot data frame
df_plot <- data.frame(
                      Investment = x_plot,
                      Actual = risk(x_plot),
                      Approximated = risk_approx(x_plot)
                     )

ggplot(df_plot, aes(x=Investment)) +
  geom_line(aes(y=Actual, color="Actual Risk")) +
  geom_line(aes(y=Approximated, color="Approximated Risk"), linetype="dashed") +
  geom_vline(xintercept = x_opt, linetype="dotted", color="darkred") +
  labs(title="Investment Risk: Actual vs Approximated",
      x="Investment Amount (x)",
      y="Risk (sqrt(x))",
      color="Legend") +
  theme_minimal()

The plot shows the accuracy of the linear approximation across a range of investments. The dotted line marks the recommended investment under the risk limit. Managers can easily see the safe region for investment using this visual.

Question 3: Inventory Management

Taylor Series Expansion

# Original cost function
C <- function(p) exp(p)

# Taylor Series expansion around p = 0 (up to second degree)
C_approx <- function(p) {
  1 + p + p^2 / 2
}

# Example: check approximation at p = 0.1
C_approx(0.1)

## [1] 1.105

C(0.1)

## [1] 1.105171

The Taylor Series approximates the exponential cost function near p = 0. Using C(p) ≈ 1 + p + 0.5*p^2 simplifies calculations and helps in quick pricing analysis without evaluating expensive exponential functions repeatedly.

Approximating Profit

# Original profit
profit_original <- function(p) {
p*(1 - p) - exp(p)
}

# Approximated profit using Taylor Series for cost
profit_approx <- function(p) {
p*(1 - p) - (1 + p + p^2 / 2)
}

# Compare profit for a range of prices
p_vals <- seq(0, 1, by=0.01)
profit_df <- data.frame(
Price = p_vals,
Original = sapply(p_vals, profit_original),
Approximation = sapply(p_vals, profit_approx)
)

# Take a look
head(profit_df, 10)

##    Price  Original Approximation
## 1   0.00 -1.000000      -1.00000
## 2   0.01 -1.000150      -1.00015
## 3   0.02 -1.000601      -1.00060
## 4   0.03 -1.001355      -1.00135
## 5   0.04 -1.002411      -1.00240
## 6   0.05 -1.003771      -1.00375
## 7   0.06 -1.005437      -1.00540
## 8   0.07 -1.007408      -1.00735
## 9   0.08 -1.009687      -1.00960
## 10  0.09 -1.012274      -1.01215

This compares original and approximated profit across prices. It helps managers see where the approximation works well, which is enough for initial pricing strategy decisions without complex computations.

Pricing Strategy

# Approximated profit: profit_approx = p*(1-p) - (1 + p + p^2/2) = p - p^2 -1 - p - p^2/2 = -3/2*p^2 -1

# Derivitive
profit_approx_prime <- function(p) -3 * p

# Solve derivative = 0
p_opt_approx <- 0
profit_max_approx <- profit_approx(p_opt_approx)

p_opt_approx

## [1] 0

profit_max_approx

## [1] -1

Using the derivative of the approximated profit, we find the price that maximizes profit in the simplified model. This method gives a fast estimate for pricing strategy. The approximate model suggests p=0 stationary-boundary solution.

Visualization

ggplot(profit_df, aes(x=Price)) +
  geom_line(aes(y=Original, color="Original Profit")) +
  geom_line(aes(y=Approximation, color="Approximated Profit"), linetype="dashed") +
  geom_vline(xintercept = p_opt_approx, linetype="dotted", color="darkred") +
  labs(title="Profit: Original vs Approximated Function",
      x="Price (p)",
      y="Profit",
      color="Function") +
  theme_minimal()

The plot shows how the approximated profit compares with the actual profit function. The dotted line indicates the estimated optimal price. This helps visually confirm that the approximation provides a reasonable guidance for pricing.

Question 4: Economic Forecasting

Maclaurin Series Expansion

# Original growth function
G_original <- function(x) log(1 + x)

# Maclaurin Series 
G_maclaurin <- function(x) x - x^2 / 2

The Maclaurin Series approximates economic growth for very small investments. It allows quick estimation of initial returns from infrastructure spending without calculating logarithms.

Approximation of Growth

# Taylor Series expansion 
x0 <- 0.2
G_taylor <- function(x) {
G0 <- G_original(x0)
G1 <- 1 / (1 + x0)
G2 <- -1 / (1 + x0)^2
G0 + G1*(x - x0) + 0.5*G2*(x - x0)^2
}

The Taylor Series gives a more accurate approximation near x0 = 0.2. It is useful for modeling growth when planning moderate infrastructure investments without recalculating logarithms each time.

# Define a range of investments
x_vals <- seq(0, 1, by = 0.01)

# Compare original vs approximated values
growth_df <- data.frame(
                        Investment = x_vals,
                        Original = sapply(x_vals, G_original),
                        Maclaurin = sapply(x_vals, G_maclaurin),
                        Taylor = sapply(x_vals, G_taylor)
)

# Show first 10 rows
head(growth_df, 10)

##    Investment    Original Maclaurin      Taylor
## 1        0.00 0.000000000   0.00000 0.001766001
## 2        0.01 0.009950331   0.00995 0.011453501
## 3        0.02 0.019802627   0.01980 0.021071557
## 4        0.03 0.029558802   0.02955 0.030620168
## 5        0.04 0.039220713   0.03920 0.040099335
## 6        0.05 0.048790164   0.04875 0.049509057
## 7        0.06 0.058268908   0.05820 0.058849335
## 8        0.07 0.067658648   0.06755 0.068120168
## 9        0.08 0.076961041   0.07680 0.077321557
## 10       0.09 0.086177696   0.08595 0.086453501

The table shows numerical differences between the actual growth and the approximations. Managers can see where each approximation is reliable and understand the limits of using simplified models.

# Plot comparison
ggplot(growth_df, aes(x = Investment)) +
  geom_line(aes(y = Original, color = "Original")) +
  geom_line(aes(y = Maclaurin, color = "Maclaurin"), linetype = "dashed") +
  geom_line(aes(y = Taylor, color = "Taylor (x0=0.2)"), linetype = "dotted") +
  labs(
      title = "Growth Function Approximations",
      x = "Investment (x)",
      y = "Growth G(x)",
      color = "Function"
      ) +
  theme_minimal()

The plot compares original growth with both approximations. The dashed line (Maclaurin) is accurate for very small investments, while the dotted line (Taylor) is more accurate near x0 = 0.2. Managers can use this visual to quickly estimate growth for different investment levels.

Policy Reccomendation

G_target <- 0.1

# Maclaurin: x - x^2/2 = G_target
a <- 0.5
b <- -1
c <- G_target

x_solution <- (-b + sqrt(b^2 - 4*a*c)) / (2*a)
x_solution

## [1] 1.894427

This shows the investment level required to reach a target growth of 0.1 using the Maclaurin approximation. While it gives a quick estimate, policymakers should remember that it is only accurate for small investments. For larger projects, using the original logarithmic function is recommended.

Problem 5: Profit, Cost, & Pricing

Question 1: Profit Maximization

# Profit function
profit_function <- function(x, y) {
  30*x - 2*x^2 - 3*x*y + 24*y - 4*y^2
}

# Partial derivatives (for reference)
partial_x <- function(x, y) 30 - 4*x - 3*y
partial_y <- function(x, y) 24 - 3*x - 8*y

# Solve linear system: 4x + 3y = 30 ; 3x + 8y = 24
A <- matrix(c(4, 3, 3, 8), nrow = 2, byrow = TRUE)
b <- c(30, 24)
solution <- solve(A, b)
x_crit <- solution[1]
y_crit <- solution[2]

# Profit at critical point
profit_crit <- profit_function(x_crit, y_crit)

# Hessian and classification (second derivative test)
d2x <- -4
d2y <- -8
d2xy <- -3
hessian_det <- d2x * d2y - d2xy^2

# Eigenvalues of Hessian (extra check)
H <- matrix(c(d2x, d2xy, d2xy, d2y), 2, 2)
hess_eigs <- eigen(H)$values

# Print a tidy numeric summary
summary_table <- data.frame(
  Quantity = c("x_crit", "y_crit", "Profit_at_crit", "Hessian_det"),
  Value = c(round(x_crit, 6), round(y_crit, 6), round(profit_crit, 6), round(hessian_det, 6))
)
kable(summary_table, caption = "Critical point summary")

Critical point summary

Quantity	Value
x_crit	7.304348
y_crit	0.260870
Profit_at_crit	112.695652
Hessian_det	23.000000

cat("\nHessian eigenvalues:", paste(round(hess_eigs, 6), collapse = ", "), "\n")

## 
## Hessian eigenvalues: -2.394449, -9.605551

# Prepare grid for plotting (z as profit)
x_vals <- seq(0, 15, length.out = 200)
y_vals <- seq(0, 10, length.out = 200)
z_matrix <- outer(x_vals, y_vals, Vectorize(function(a, b) profit_function(a, b)))

The profit function has one critical point, located at (x = 7.304, y = 0.261). The Hessian test shows this point is a local maximum, meaning the company earns the highest profit in that region when producing these quantities. At this point, the profit is PI = 112.696. There are no local minima or saddle points, so this is the only meaningful optimal production level for the firm.

Discussion: The results show that the company maximizes profit by producing about 7.3 units of product A and 0.26 units of product B. This production mix gives the highest possible profit of roughly 113. There are no saddle points, so small changes around this level clearly reduce profit. The company should focus on maintaining production near this optimal point to avoid losses. Any deviations should be carefully monitored, but the risk from nearby points is low since the maximum is well-defined.

Question 2: Pricing Strategy

# Demand functions
demand_x <- function(x, y) 120 - 15*x + 10*y
demand_y <- function(x, y) 80 + 5*x - 20*y

# Revenue function
revenue_function <- function(x, y) {
  x * demand_x(x, y) + y * demand_y(x, y)
}

# Partial derivatives
revenue_partial_x <- function(x, y) 120 - 30*x + 15*y
revenue_partial_y <- function(x, y) 80 + 15*x - 40*y

# Solve for critical point
A <- matrix(c(-30, 15, 15, -40), nrow = 2, byrow = TRUE)
b <- c(-120, -80)
solution <- solve(A, b)

x_crit <- solution[1]
y_crit <- solution[2]

# Revenue at optimal prices
revenue_crit <- revenue_function(x_crit, y_crit)

# Demands at optimal prices
dx_opt <- demand_x(x_crit, y_crit)
dy_opt <- demand_y(x_crit, y_crit)

# Hessian determinant to classify critical point
d2x <- -30
d2y <- -40
d2xy <- 15
hessian_det <- d2x * d2y - d2xy^2

## Optimal Brand X Price: $ 6.15

## Optimal Brand Y Price: $ 4.31

## Maximum Revenue: $ 541.54

## Brand X Demand: 70.77 units

## Brand Y Demand: 24.62 units

## Hessian Determinant: 975

Discussion: This analysis finds the optimal pricing strategy for two competing brands based on how price changes affect demand. The revenue function has one critical point, and the Hessian test confirms it is a local maximum, meaning it produces the highest possible revenue. The store maximizes revenue by setting Brand X at 6.15 and Brand Y at 4.31, resulting in a total revenue of roughly 542. Because the Hessian determinant is positive and both second derivatives are negative, there are no saddle points, giving the manager a clear and stable pricing recommendation.

Question 3: Cost Minimization

# Cost function
cost_function <- function(x, y) {
  18*x^2 + 110*y^2 + 12*x + 18*y + 1500
}

# Cost function with constraint substituted
cost_constrained <- function(x) {
  y <- 200 - x
  cost_function(x, y)
}

# Optimal solution (first derivative = 0)
x_opt <- 44006 / 256
y_opt <- 200 - x_opt

# Minimum cost
min_cost <- cost_function(x_opt, y_opt)

## Optimal New York production: 171.9 units

## Optimal Chicago production: 28.1 units

## Minimum total cost: $ 622818.7

Discussion: This analysis determines how a company should split production between its New York and Chicago plants to minimize weekly costs. The optimal allocation is 171.90 units in New York and 28.10 units in Chicago, resulting in a minimized total cost of approximately $622,819. Although New York has slightly higher linear costs, Chicago’s much larger quadratic cost term makes high production there disproportionately expensive. Concentrating production in New York is the most economical strategy, though real-world considerations (such as labor constraints, factory capacity, or shipping logistics) may require adjustments.

Question 4: Marketing Mix

# Effectiveness function
effectiveness_function <- function(x, y) {
  500*x + 700*y - 5*x^2 - 10*x*y - 8*y^2
}

# Partial derivatives (for reference)
effective_partial_x <- function(x, y) 500 - 10*x - 10*y
effective_partial_y <- function(x, y) 700 - 10*x - 16*y

# Solve analytically
y_opt <- 200 / 6
x_opt <- 50 - y_opt

# Evaluate maximum effectiveness
max_effectiveness <- effectiveness_function(x_opt, y_opt)

# Hessian check
d2x <- -10
d2y <- -16
d2xy <- -10
hessian_det <- d2x*d2y - d2xy^2

## Optimal online ad spending: 16.67 thousand dollars

## Optimal TV ad spending: 33.33 thousand dollars

## Maximum customer reach: 15833.33

## Hessian determinant: 60

Discussion: The results suggest putting more money into the channel that gives the biggest customer reach, roughly 2/3 to TV and 1/3 to online, as a starting point. Treat this allocation as a pilot and measure performance before scaling up. Track key metrics like reach, clicks, and conversions weekly, and adjust quickly if one channel underperforms. If a saddle point occurs, it means small changes don’t clearly improve results, so rely on short experiments to test nearby budget splits. Add practical constraints and gather more data if needed, then adjust the allocation based on real-world performance.