Your grandparents have an annuity. The value of the annuity increases each month by an automatic deposit of 1% interest on the 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.
\[ a_{n+1} = 1.01a_n - 1000 \]
Will the annuity run out of money? When?
# N0
amount = 50000
# Counter for the loop
i <- 0
# loop through periods in dynamical model while balance > 0
while (amount > 0){
amount = 1.01*amount - 1000
print(amount)
i <- i + 1
print(i)
if (amount < 1000) break # break out of loop if balance is less than 1000
}
## [1] 49500
## [1] 1
## [1] 48995
## [1] 2
## [1] 48484.95
## [1] 3
## [1] 47969.8
## [1] 4
## [1] 47449.5
## [1] 5
## [1] 46923.99
## [1] 6
## [1] 46393.23
## [1] 7
## [1] 45857.16
## [1] 8
## [1] 45315.74
## [1] 9
## [1] 44768.89
## [1] 10
## [1] 44216.58
## [1] 11
## [1] 43658.75
## [1] 12
## [1] 43095.34
## [1] 13
## [1] 42526.29
## [1] 14
## [1] 41951.55
## [1] 15
## [1] 41371.07
## [1] 16
## [1] 40784.78
## [1] 17
## [1] 40192.63
## [1] 18
## [1] 39594.55
## [1] 19
## [1] 38990.5
## [1] 20
## [1] 38380.4
## [1] 21
## [1] 37764.21
## [1] 22
## [1] 37141.85
## [1] 23
## [1] 36513.27
## [1] 24
## [1] 35878.4
## [1] 25
## [1] 35237.18
## [1] 26
## [1] 34589.56
## [1] 27
## [1] 33935.45
## [1] 28
## [1] 33274.81
## [1] 29
## [1] 32607.55
## [1] 30
## [1] 31933.63
## [1] 31
## [1] 31252.97
## [1] 32
## [1] 30565.5
## [1] 33
## [1] 29871.15
## [1] 34
## [1] 29169.86
## [1] 35
## [1] 28461.56
## [1] 36
## [1] 27746.18
## [1] 37
## [1] 27023.64
## [1] 38
## [1] 26293.87
## [1] 39
## [1] 25556.81
## [1] 40
## [1] 24812.38
## [1] 41
## [1] 24060.51
## [1] 42
## [1] 23301.11
## [1] 43
## [1] 22534.12
## [1] 44
## [1] 21759.46
## [1] 45
## [1] 20977.06
## [1] 46
## [1] 20186.83
## [1] 47
## [1] 19388.7
## [1] 48
## [1] 18582.58
## [1] 49
## [1] 17768.41
## [1] 50
## [1] 16946.09
## [1] 51
## [1] 16115.55
## [1] 52
## [1] 15276.71
## [1] 53
## [1] 14429.48
## [1] 54
## [1] 13573.77
## [1] 55
## [1] 12709.51
## [1] 56
## [1] 11836.6
## [1] 57
## [1] 10954.97
## [1] 58
## [1] 10064.52
## [1] 59
## [1] 9165.165
## [1] 60
## [1] 8256.817
## [1] 61
## [1] 7339.385
## [1] 62
## [1] 6412.779
## [1] 63
## [1] 5476.907
## [1] 64
## [1] 4531.676
## [1] 65
## [1] 3576.992
## [1] 66
## [1] 2612.762
## [1] 67
## [1] 1638.89
## [1] 68
## [1] 655.2788
## [1] 69
Starting with month 1, the annuity will end up with $655.2788 in month 69.
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\times 5 = 30mph\)) requires a stopping distances of \(a_6 = 47ft\)
library(ggplot2)
# DaTA frame for the data in table
bDistancedf <- data.frame(n=1:16, an=c(3,6,11,21,32,47,65,87,112,140,171,204,241,282,325,376))
# Calc Delta N
bDistancedf$deltaAn = c(NA,tail(bDistancedf$an,-1)-head(bDistancedf$an,-1))
#Plot
ggplot(bDistancedf, aes(x=n, y=deltaAn)) + geom_point()
## Warning: Removed 1 rows containing missing values (geom_point).
The graph does reasonably approximate a linear relationship.
Since this is linear, we’re trying to find a line. We can also assume that if \(n = 0\), the stopping distance would be \(0\), and therefore the line passes through the origin.
To approximate the slope of this line, lets average the slopes from the origin to each of the points. In other words, for each point, lets calculate \(\frac{\Delta a_n}{n}\)
# Calc Slope y/x
bDistancedf$slope <- bDistancedf$deltaAn / bDistancedf$n
slope <- mean(bDistancedf$slope, na.rm=TRUE)
We could easily find the least squares solution to this problem, forcing the intercept to be 0, but the textbook was solving these problems in an approximate way by calculating the slope
altSlope <- summary(lm(deltaAn ~ 0 + n, data=bDistancedf))$coefficients["n","Estimate"]
Lets plot the errors in the predicted values vs n
bDistancedf$estimate <- slope*bDistancedf$n
bDistancedf$error <- bDistancedf$estimate - bDistancedf$deltaAn
ggplot(bDistancedf, aes(x=n, y=error)) + geom_point()
## Warning: Removed 1 rows containing missing values (geom_point).
This graph suggests that there is a pattern to the errors. This model may not be appropriate. If we were just fitting \(\Delta a_n\) to \(n\) using a linear model including an intercept, we could perhaps end up with a better fit.
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 (example 3, section 1.2) in that the number of people hearing the rumor each day is proportional to the product of the number 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 \(k\) is a parameter that depends on how fast the rumor spreads and \(n\) is the number of days. Assume \(k = 0.001\) and further assume that four people initially have heard the rumor. How soon will all 1000 employees have heard the rumor?
rumordf <- data.frame(n=0, r=4, roundr=4)
k <- 0.001
for(i in 1:13){
rumordf <- rbind(rumordf, c(i, rumordf$r[i] + k*rumordf$r[i]*(1000 - rumordf$r[i]),
round(rumordf$roundr[i] + k*rumordf$roundr[i]*(1000 -
rumordf$roundr[i]))))
}
rumordf
## n r roundr
## 1 0 4.00000 4
## 2 1 7.98400 8
## 3 2 15.90426 16
## 4 3 31.55557 32
## 5 4 62.11538 63
## 6 5 120.37244 122
## 7 6 226.25535 229
## 8 7 401.31922 406
## 9 8 641.58132 647
## 10 9 871.53605 875
## 11 10 983.49701 984
## 12 11 999.72765 1000
## 13 12 999.99993 1000
## 14 13 1000.00000 1000
I created two columns here, one that rounds the number of people who have heard the rumor to the nearest integer. Given the initial condition of four people who have heard the rumor, this doesn’t make too big a difference. After 12 days, the rumor will have spread throughout the office.
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. Over time, there is an interaction between price and supply. The economist has proposed the following model,p 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.
\[ P_{n+1} = P_n - 0.1(Q_n - 500) \] \[ Q_{n+1} = Q_n + 0.2(P_n - 100) \]
See explanation below.
Case | Price | Quantity |
---|---|---|
A | 100 | 500 |
B | 200 | 500 |
C | 100 | 600 |
D | 100 | 400 |
library(reshape2)
iterations <- 200
pinit <- 100
qinit <- 500
econdf <- data.frame(n = 0, p = pinit, q= qinit)
for(i in 1:iterations){
econdf <- rbind(econdf,
c(i, tail(econdf$p,1) - 0.1*(tail(econdf$q,1) - 500),
tail(econdf$q,1) + 0.2*(tail(econdf$p,1) - 100)))
}
ggplot(melt(econdf, id="n"), aes(x=n, y=value, color=variable)) + geom_line()
#################################################################################
pinit <- 200
qinit <- 500
econdf <- data.frame(n = 0, p = pinit, q= qinit)
for(i in 1:iterations){
econdf <- rbind(econdf,
c(i, tail(econdf$p,1) - 0.1*(tail(econdf$q,1) - 500),
tail(econdf$q,1) + 0.2*(tail(econdf$p,1) - 100)))
}
ggplot(melt(econdf, id="n"), aes(x=n, y=value, color=variable)) + geom_line()
#################################################################################
pinit <- 100
qinit <- 600
econdf <- data.frame(n = 0, p = pinit, q= qinit)
for(i in 1:iterations){
econdf <- rbind(econdf,
c(i, tail(econdf$p,1) - 0.1*(tail(econdf$q,1) - 500),
tail(econdf$q,1) + 0.2*(tail(econdf$p,1) - 100)))
}
ggplot(melt(econdf, id="n"), aes(x=n, y=value, color=variable)) + geom_line()
#################################################################################
pinit <- 100
qinit <- 400
econdf <- data.frame(n = 0, p = pinit, q= qinit)
for(i in 1:iterations){
econdf <- rbind(econdf,
c(i, tail(econdf$p,1) - 0.1*(tail(econdf$q,1) - 500),
tail(econdf$q,1) + 0.2*(tail(econdf$p,1) - 100)))
}
ggplot(melt(econdf, id="n"), aes(x=n, y=value, color=variable)) + geom_line()
Seeing the graphs made me realize what was going on intuitively. At first I assumed this system would result in an equilibrium price and quantity. While this is possible given the right initial conditions (a price of 100 and a quantity of 500), any other values will cause the system to oscillate without bound around the equilibrium.
The difference in sign between 0.1 and 0.2 in the two equations is what causes the oscillation. If these were both positive or both negative, they would increase or decrease without bound most of the time (this depends on what the initial price and quantity are).
The one equilibrium case I found was if the equations are the same, the second term is being subtracted from the first, and the initial price and quantity are the same:
pinit <- 20
qinit <- 20
econdf <- data.frame(n = 0, p = pinit, q= qinit)
for(i in 1:iterations){
econdf <- rbind(econdf,
c(i, tail(econdf$p,1) - 0.1*(tail(econdf$q,1) - 500),
tail(econdf$q,1) - 0.1*(tail(econdf$p,1) - 500)))
}
ggplot(melt(econdf, id="n"), aes(x=n, y=value, color=variable)) + geom_line()
Any differences from this specific set of cases will once again cause divergence. For example, in the above case if the initial price is not equal to the initial quantity:
pinit <- 21
qinit <- 20
econdf <- data.frame(n = 0, p = pinit, q= qinit)
for(i in 1:iterations){
econdf <- rbind(econdf,
c(i, tail(econdf$p,1) - 0.1*(tail(econdf$q,1) - 500),
tail(econdf$q,1) - 0.1*(tail(econdf$p,1) - 500)))
}
ggplot(melt(econdf, id="n"), aes(x=n, y=value, color=variable)) + geom_line()
This is probably not a good model for the interaction between price and quantity in normal markets. However, I could imagine some economic processes that might function in this way. This sort of oscillating behavior reminded me of the “fool in the shower” metaphor for monetary policy, where central banks with imperfect information react to the economy in a delayed fashion, raising and lowering interest rates in a way that’s out of sync with the rest of the economy.