%% Lab 5 -- Boom Bust Butterfly
%% A model of the market
% As described in the pre-lab, we will be modeling the number of producers in a given industry X(t) where t is the number of financial quarters since the first quarter of 1980.
% We use a response curve R(X) to relate the number of insutry players in the next quarter to the number active in the current quarter by the following formula
% X(t + 1) = R(X(t)) \cdot X(t).
% See the pre-lab for an example of how this relationship defines a numerical sequence.
%
%% Exercise 1
% Typically the rate of growth function R(X) exhibits the following features:
% R(X) is a decreasing function of X.
% R(X) is never less than zero.
% In the context of our model (market economy) described in the pre-lab, explain why, intuitively, these two features should be expected. Also, explain briefly why the solution X_0 to R(X) = 1 is the market equilibrium.
%
%% Volatility
% In this lab we will look at the same market but with different degrees of volatility. We will denote the volatility index by v. A market with higher volatility will show more extreme responses:
% When the number of producers is below market equlibrium: the higher the volatility, the more new producers will enter the market.
% When the number of producers is above market equilibrium: the higher the volatility, the faster existing producers will exit the market.
% Our model will use the following response curve, given as a function both of the current number of producers X and the volatility index v.  (The precise function is not that important; it is the qualitative features discussed in the above paragraphs that counts. We give an explicit function here so we can run simulations later.)
%   R(X,v) = \sqrt{1 + v^2(x-1)^2} (1 - \tanh(v(x-1)))
% Let's plot it and check that it captures all the qualitative features described previously.

% Define our response curve
R = @(X,v) sqrt(1 + v^2.*(X-1).^2) .* (1 - tanh(v .* (X-1)))

% Plot it for a few different values of volatility v
Xvalues = linspace(0,5,200);
plot(Xvalues,R(Xvalues,1), 'r') % Red is v = 1
hold on;
plot(Xvalues, R(Xvalues,2), 'g') % Green is v = 2
plot(Xvalues,R(Xvalues,0.5), 'b') % Blue is v = 0.5
set(gca, 'FontSize', 16)
xlabel({'X, number of producers', '(measured in hundreds)'});
ylabel('R, response rate');
hold off;

% We see indeed:
% All three response curves are decreasing.
% They model the same market: they have the same market equilibrium. (Which we have chosen to be at X = 1 (corresponding to 100 producers)).
% Higher volatility gives more extreme responses to deviations from equilibrium
%
%% Exercise 2
% When X(0) = 0.4, calculate the percent change from X(0) to X(1) for the volatility levels v = 0.5, 1, and 2.
% Do the same as above, but with starting value X(0) = 2.
% You may use MATLAB to aid in completing this part.

X0 = 0.5
X1 = R(X0,1)* X0
% Above computes X(1) when X(0) = 0.5 and volatility v = 1. 
% How do you use this to compute the percent change from X0 to X1?

%
%% Simulations
% Now let's see some simulations.
%% Case 1: Low volatility
% The code below plots the trajectory X(t) for three randomly generated different starting values X(0), when the volatility v is low (here low volatility is when 0 < v < 1).

v = 0.4;   % Set volatility to a small number

MaxNumQuarters = 20; % Carry out the simulation for how many quarters

Xa = zeros(MaxNumQuarters,1);  % Pre-allocate empty arrays to speed
Xb = zeros(MaxNumQuarters,1);  %  up computation
Xc = zeros(MaxNumQuarters,1);

Xa(1) = rand; % Generate three different initial values randomly
Xb(1) = rand; %    between 0 and 1 
Xc(1) = rand;

for t = 1:MaxNumQuarters-1   % Build the time-series corresponding to the values
                % using the formula X(t+1) = R(X(t), v) X(t)
    Xa(t+1) = R(Xa(t), v) * Xa(t);
    Xb(t+1) = R(Xb(t), v) * Xb(t);
    Xc(t+1) = R(Xc(t), v) * Xc(t);
end

plot(Xa, 'r.-')
hold on;
plot(Xb, 'g.-')
plot(Xc, 'b.-')
hold off;

% You will see that in all the cases shown, the number of companies in the market steadily increases and converges toward the market equilibrium.
%
%% Exercise 3
% In the low volatility setting, does the convergence toward market equilibrium happen faster for less volatile markets (e.g. v = 0.2) or more volatile markets (e.g. v = 0.8)? How does this correspond to the notion of volatility (as described above in this lab)?
%
%% Case 2: Medium volatility
% Something different happens when the market volatility becomes a medium value. For our purposes this is when 1 < v < 2. Run the simulation below to see what happens.

v = 1.4;   % Set volatility to a larger number

MaxNumQuarters = 20; % Carry out the simulation for how many quarters

Xa = zeros(MaxNumQuarters,1);  % Pre-allocate empty arrays to speed
Xb = zeros(MaxNumQuarters,1);  %  up computation
Xc = zeros(MaxNumQuarters,1);

Xa(1) = rand; % Generate three different initial values randomly
Xb(1) = rand; %    between 0 and 1 
Xc(1) = rand;

for t = 1:MaxNumQuarters-1   % Build the time-series corresponding to the values
                % using the formula X(t+1) = R(X(t), v) X(t)
    Xa(t+1) = R(Xa(t), v) * Xa(t);
    Xb(t+1) = R(Xb(t), v) * Xb(t);
    Xc(t+1) = R(Xc(t), v) * Xc(t);
end

plot(Xa, 'r.-')
hold on;
plot(Xb, 'g.-')
plot(Xc, 'b.-')
hold off;

% What you will see is that, in all cases, the number of companies in the market is no longer always increasing. The number of companies sometimes overshoots the market equilibrium. But in all cases eventually the number of companies converges to market equilibrium.
%
%% Exercise 4
% In this medium volatility setting, does the convergence toward market equilibrium happen faster for less volatile markets (e.g. v = 1.2) or more volatile markets (e.g. v = 1.8)? How does this correspond to the notion of volatility (as described above in this lab)?
%
%% Case 3: High volatility
% When the volatility exceeds 2, something very different happens.
% First, when v is between 2 and approximately 2.4, this happens.

v = 2.3;   % Set volatility to a larger number

MaxNumQuarters = 30; % Carry out the simulation for how many quarters

Xa = zeros(MaxNumQuarters,1);  % Pre-allocate empty arrays to speed
Xb = zeros(MaxNumQuarters,1);  %  up computation
Xc = zeros(MaxNumQuarters,1);

Xa(1) = rand; % Generate three different initial values randomly
Xb(1) = rand; %    between 0 and 1 
Xc(1) = rand;

for t = 1:MaxNumQuarters-1   % Build the time-series corresponding to the values
                % using the formula X(t+1) = R(X(t), v) X(t)
    Xa(t+1) = R(Xa(t), v) * Xa(t);
    Xb(t+1) = R(Xb(t), v) * Xb(t);
    Xc(t+1) = R(Xc(t), v) * Xc(t);
end

plot(Xa, 'r.-')
hold on;
plot(Xb, 'g.-')
plot(Xc, 'b.-')
hold off;

% When v is in this range (between 2 and approximately 2.4), we see that
% The number of companies no longer converge to market equilibrium!
% Instead we have the new phenomenon where
% The number of companies settles into a boom-bust cycle.
% Irregardless of the starting number of companies, eventually the market settles into a regular pattern where a quarter with large number of producers is followed by a quarter with a small number of producers.
% What if we increase v even futher?
% When v is around 2.5 the market settles into a four quarter cycle, with
% First a very good quarter (lots of producers)
% Next a very bad quarter (very few producers)
% Next a somewhat good quarter (many producers, but not as many as the first quarter)
% Next a somewhat bad quarter (few producers, but not as few as the second quarter)
% repeat.

v = 2.5;   % Set volatility to a larger number

MaxNumQuarters = 30; % Carry out the simulation for how many quarters

Xa = zeros(MaxNumQuarters,1);  % Pre-allocate empty arrays to speed
Xb = zeros(MaxNumQuarters,1);  %  up computation
Xc = zeros(MaxNumQuarters,1);

Xa(1) = rand; % Generate three different initial values randomly
Xb(1) = rand; %    between 0 and 1 
Xc(1) = rand;

for t = 1:MaxNumQuarters-1   % Build the time-series corresponding to the values
                % using the formula X(t+1) = R(X(t), v) X(t)
    Xa(t+1) = R(Xa(t), v) * Xa(t);
    Xb(t+1) = R(Xb(t), v) * Xb(t);
    Xc(t+1) = R(Xc(t), v) * Xc(t);
end

plot(Xa, 'r.-')
hold on;
plot(Xb, 'g.-')
plot(Xc, 'b.-')
hold off;

% And after you increase v above 2.6, the pattern becomes harder and harder to discern. You start picking up solutions with all different sorts of periods.

v = 2.7;   % Set volatility to a larger number

MaxNumQuarters = 30; % Carry out the simulation for how many quarters

% Just one graph now, to make it clearer what's happening
Xa = zeros(MaxNumQuarters,1);  % Pre-allocate empty array to speed up
Xb = zeros(MaxNumQuarters,1);

Xa(1) = rand; % Generate random starting value
Xb(1) = rand; 

for t = 1:MaxNumQuarters-1   % Build the time-series corresponding to the values
                % using the formula X(t+1) = R(X(t), v) X(t)
    Xa(t+1) = R(Xa(t), v) * Xa(t);
    Xb(t+1) = R(Xb(t), v) * Xb(t);
end

plot(Xa, 'r.-')
hold on;
plot(Xb, 'b.-')
hold off;

% This is the FIRST reason why it is hard to predict the stock market: the patterns are just not as simple as booms and busts. When the volatility is high, the same market conditions can support cycles of differing periods.
%
%% Chaos
% There is a second difference between highly volatile and less volatile markets: the famed butterfly effect. You can see it in action below.
% The code below allows you to simulate four different settings:
% The low volatility setting when v = 0.4
% The medium volatility setting when v = 1.8
% The high volatitliy boom-bust scenario when v = 2.2
% The high volatility setting when v = 2.8
% The simulation starts out with two initial values that are close to each other: they differ by 5%.

v = 2.8;   % Set volatility 

MaxNumQuarters = 30; % Carry out the simulation for how many quarters

% Just one graph now, to make it clearer what's happening
Xa = zeros(MaxNumQuarters,1);  % Pre-allocate empty array to speed up
Xb = zeros(MaxNumQuarters,1);

Xa(1) = rand; % Generate random starting value
Xb(1) = Xa(1)*1.05;   % The second starting value differs from the first by 5%

for t = 1:MaxNumQuarters-1   % Build the time-series corresponding to the values
                % using the formula X(t+1) = R(X(t), v) X(t)
    Xa(t+1) = R(Xa(t), v) * Xa(t);
    Xb(t+1) = R(Xb(t), v) * Xb(t);
end

plot(Xa, 'r.-')
hold on;
plot(Xb, 'b.-')
hold off;

%
%% Exercise 5
% Run the simulation in the "Chaos" section a few times, with different values of volatility. What do you notice that is different qualitatively between the v = 2.8 case versus the other three cases? What does this say about attempts to forecast stock prices, in a very highly volatile market, based on incomplete information?
%
%% Some further notes
% In our lab we chose a very specific response curve. The type of phenomea observed in the lab is extremely prevalent and just depends on a few qualitative features of the response curve. Generally speaking, as long as the response curve exhibits a "hump", the boom-bust cycle would appear when the volatility is high.
% Something interesting happens when v becomes 3 or larger. When v = 3 this marks the appearance of the "cycle with period 3". Mathematical research has shown that when the volatility is high enough so cycles of period 3 shows up, depending on the initial condition, cycles of all possible periods will also show up. Meaning that once you get to v = 3 and larger, the system becomes as unpredictable as possible.
%% References
% The cyclical booms and busts, the dependence on the volatility, and the appearance of chaos, are very common features in the biological and social sciences in addition to economics. If you want to learn a bit more about this, here are some approachable references.
% The Nature article Simple Mathematical Models with very Complicated Dynamics by Robert May (Nature 261, 459–467) surveys the basic mathematical ideas as well as some of the applications.
% W.E. Ricker's Stock and Recruitment (Journal of the Fisheries Research Board of Canada, 1954, 11(5): 559-623) examines how environmental pressure and human interference can combine to lead to cyclical booms and busts in fish population.
% William J. Baumol and Jess Benhabib's Chaos: Significance, Mechanism, and Economic Applications (The Journal of Economic Perspectives 3, 77-105) surveys the many ways chaotic behavior can arise in simple economic models.
% The mathematical theory describing how the system transitions from a relatively mild one at low volatility to a completely unpredictable one at high volatility is called chaos theory. A very good account of chaos theory, emphasizing on the history of its discovery and scientific implications (with only a little dose of mathematics), can be found in
% James Gleick, Chaos: Making a New Science, Viking Books. (eBook version: Open Road Media)