DATA 609: Week 1 Homework
Keith Folsom
February 5, 2017
suppressWarnings(suppressMessages(library(ggplot2)))
suppressWarnings(suppressMessages(library(dplyr)))Page 8: #10
Your parents have an annuity. The value of the annuity increases each month by an automatci deposit of 1% interest on th previous month’s balance. Your grandparents withdraw $1000 at the beginning of each month for living expenses. Currently they have $50,000 in the annuity. Model the annuity with a dynamical system. Will the annuity run out of money?
\[a_{n + 1} = a_n * 0.01 - 1000, a_0 = 50000\]
annuity_model <- function(amt) {
return (amt * 0.01 - 1000)
}
annuity_amt <- 50000
n <- 1
while (TRUE) {
annuity_amt[n+1] <- annuity_amt[n] + annuity_model(annuity_amt[n])
if (annuity_amt[n+1] < 0) break
n <- n + 1
}
df <- data.frame(month = 1:length(annuity_amt), amount = annuity_amt)
ggplot(df, aes(x=month, y=amount)) + geom_line() + scale_x_continuous(limits = c(0, 80)) + ggtitle("Annuity Amount Over Time")
The annuity will run out of money after month 70.
## month amount
## 67 67 3576.9923
## 68 68 2612.7623
## 69 69 1638.8899
## 70 70 655.2788
## 71 71 -338.1684Page 17: #9
The data in the accompanying table show the speed \(n\) (in increments of 5 mph) of an automobile and the associated distance \(a_n\) in feet required to stop it once the brakes are applied. For instance, \(n = 6\) ) representing 6 x 5 = 30 mph) requires a stopping distance of \(a_6\) = 47 ft.
n <- seq(1:16)
an <- c(3, 6, 11, 21, 32, 47, 65, 87, 112, 140, 171, 204, 241, 282, 325, 376)
df <- data.frame(n, an)
# --------------------
# calculate delta n
# --------------------
# for (i in 1:nrow(df)) {
# df$delta[i] <- df$an[i + 1] - df$an[i]
# }
df$delta_n <- c(diff(an), NA)
qplot(n, delta_n, data = na.omit(df), xlab = "n", ylab = "Delta n") 
The resulting graph does reasonably approximate a linear relationship.
calculate the slope:
slope <- mean(df$delta_n / df$n, na.rm = T)
slope## [1] 3.051394df$predicted <- slope * df$n
rbind(data.frame(value = 'Actual', x = df$n, y = df$an),
data.frame(value = 'Predicted', x = df$n, y = df$predicted)) %>%
ggplot(aes(x=x, y=y, color=value)) + geom_line()
ggplot(df, aes(x=n, y=(an - predicted))) + geom_point() + labs(y="Error Amount")
The model above does not result in an accurate fit especially as \(n\) increases. We see the error amount increase significantly as \(n\) increases.
Modeling \(a_n\) against n does not result in a better model either:
slope <- mean(df$an / df$n, na.rm = T)
slope## [1] 12.00929df$predicted <- slope * df$n
rbind(data.frame(value = 'Actual', x = df$n, y = df$an),
data.frame(value = 'Predicted', x = df$n, y = df$predicted)) %>%
ggplot(aes(x=x, y=y, color=value)) + geom_line()
ggplot(df, aes(x=n, y=(an - predicted))) + geom_point() + labs(y="Error Amount")
Neither model seems appropriate for this data. Modeling this as a linear model is not sufficient to create an accurate model. We see that as \(n\) increases the error in both models increases.
Page 34: #13
Consider the spreading of a rumor through a company of 1000 employees, all working in the same building. We assume that the spreading of a rumor is similar to the spreading of a contagious disease in the the number of people hearing the rumor each day is proportional to the product of the number of who have heard the rumor previously and the number who have not heard the rumor. This is given by
\[r_{n+1} = r_n + kr_n(1000 - n)\]
where is a parameter that depends on how fast the rumor spreads and \(n\) is the number of days. Assum \(k\) = 0.001 and further assume that four people initially have heard the rumor. How soon will all 1000 employees have heard the rumor?
Note. For this problem, I am assuming that the equation provided has a typo and should be:
\[r_{n+1} = r_n + kr_n(1000 - r_n)\]
k <- 0.001 # how fast the rumor spreads
n <- 1 # number of days, starting at 1
rumor_start <- 4 # number of people to initially hear the rumor
total_population <- 1000
total_rumor <- rumor_start
while (total_rumor[n] < total_population) {
total_rumor[n + 1] <- total_rumor[n] + (k*(total_rumor[n] ) * (1000 - total_rumor[n] ))
n <- n + 1
}
result <- data.frame(days = seq(1:n), number_heard_rumor = total_rumor)
result## days number_heard_rumor
## 1 1 4.00000
## 2 2 7.98400
## 3 3 15.90426
## 4 4 31.55557
## 5 5 62.11538
## 6 6 120.37244
## 7 7 226.25535
## 8 8 401.31922
## 9 9 641.58132
## 10 10 871.53605
## 11 11 983.49701
## 12 12 999.72765
## 13 13 999.99993
## 14 14 1000.00000
## 15 15 1000.00000Based on this assumption, it looks like by day 12 all 1,000 employees will have heard the rumor.
Page 55: #6
An economist is interested in the variation of the price of a single product. It is observed that a high price for the product in the market attracts more suppliers. However, increasing the quantity of the product supplied tends to drive the price down. Overtime, there is an interaction between price and supply. The economist has proposed the following model, where \(P_n\) represents the price of the product at year \(n\), and \(Q_n\) represents the quantity. Find the equilibrium values for this system.
This model does appear to be intuitive, especially in light of the relationship between price and supply. The constant of 500 in the first equation serves to keep the price positive as long as the quantity is not greater than 500. If the price exceeds 500, the price will begin decline.
Similarly, in the second equation, as price exceeds 100 the quantity will increase presumably as demand weakens.
The constant -0.1 is significant because it drops the price as the quantity increases past 500. The constant 0.2 does the reverse for quantity – as the price increases, the quantity increases.
The equilibrium points will be 500 and 100 for the quantity and price.
price_model <- function(p, q) {
return (p - (0.1 * (q - 500)))
}
quantity_model <- function(p, q) {
return (q + (0.2 * (p - 100)))
}
# ====================================================
# create a generic function to test Cases A - B
# ====================================================
test <- function(p_start, q_start, years) {
p <- p_start
q <- q_start
for (i in 1:years) {
p[i + 1] <- price_model(p[i], q[i])
q[i + 1] <- quantity_model(p[i], q[i])
}
return (data.frame(year = seq(1:(years +1)), price = p, quantity = q))
}Case A: Price = 100; Quantity = 500
caseA <- test(100, 500, 50)
ggplot(caseA, aes(year)) +
geom_line(aes(y = price, colour = "price")) +
geom_line(aes(y = quantity, colour = "quantity")) + ggtitle("Case A")
Long-term behavior shows that the quantity remains flat at 500 and the price remains flat at 100.
Case B: Price = 200; Quantity = 500
caseB <- test(200, 500, 20)
ggplot(caseB, aes(year)) +
geom_line(aes(y = price, colour = "price")) +
geom_line(aes(y = quantity, colour = "quantity")) + ggtitle("Case B")
Long-term behavior shows that the initial price is higher than the threshold value of 100 so the demand is low. As the supply or quantity increases as a consequence, the price falls dramatically.
Case C: Price = 100; Quantity = 600
caseC <- test(100, 600, 30)
ggplot(caseC, aes(year)) +
geom_line(aes(y = price, colour = "price")) +
geom_line(aes(y = quantity, colour = "quantity")) + ggtitle("Case C")
Long-term behavior shows the price starting at 100 and falling since the quantity is above 500. However, equilibrium is not reached so the long-term behavior shows the cyclic and associated rise and fall of quantity vs. price.
Case D: Price = 100; Quantity = 400
caseD <- test(100, 400, 50)
ggplot(caseD, aes(year)) +
geom_line(aes(y = price, colour = "price")) +
geom_line(aes(y = quantity, colour = "quantity")) + ggtitle("Case D")
Long-term we see an price increase leading to higher quantity and ultimately a dramatic price drop. This system does not equalize long term.