SHOOTING ALGORITHM

Part a) Fully define the competitive equilibrium

A competitive equilibrium consists of A) prices , B) allocations for the household , C) allocations for the firm  and D) firm and household profits , such that:
1) Given prices A) , the allocation of the household B)  solves the household problem 
2) Given prices A) , the allocation of the firm C)  solves the firm problem  where, 
3) Markets clear: 
				(Goods Market)
				(Labor Market)
				(Capital Market)
4) Firm profits equal household profits:

Part b. Find the first order conditions. Be sure to show all working.

The set up for the problem is as follows:
The representative household solves:
						
subject to:
						
The representative firm solves:
						
where, 
Prices adjust so that markets clear:
						(Goods Market)
						(Labor Market)
						(Capital Market)
To find the FOC for the household we will substitute the budget constraint into the utility function. From the capital accumulation constraint rearrange to get (Subscript[I, t])= (Subscript[K, t+1])-(1-δ)(Subscript[K, t]). Substitute (Subscript[I, t]) into   to get 
(Subscript[C, t]) + (Subscript[K, t+1])-(1-δ)(Subscript[K, t])= (Subscript[w, t])(Subscript[L, t]) + (Subscript[r, t])(Subscript[K, t]) + (Subscript[π, t])
(Subscript[C, t]) = (Subscript[w, t])(Subscript[L, t]) + (Subscript[r, t])(Subscript[K, t]) + (Subscript[π, t]) - (Subscript[K, t+1])+(1-δ)(Subscript[K, t])
(Subscript[C, t]) = (Subscript[w, t])(Subscript[L, t]) +(1+ (Subscript[r, t])-δ)(Subscript[K, t]) + (Subscript[π, t]) - (Subscript[K, t+1])
(Subscript[L, t]) = 1 so,
(Subscript[C, t]) = (Subscript[w, t]) +(1+ (Subscript[r, t])-δ)(Subscript[K, t]) + (Subscript[π, t]) - (Subscript[K, t+1])
Now substitute (Subscript[C, t]) into the utility function to get:
(∑_(t=0)^∞)(β^t)log((Subscript[w, t])+(1+(Subscript[r, t])-δ)(Subscript[K, t])+(Subscript[π, t])-(Subscript[K, t+1]))
Now differentiate with respect to (Subscript[K, t+1]). Notice that only terms (Subscript[C, t]) and (Subscript[C, t+1]) contain the term (Subscript[K, t+1]) so from the entire sum we can focus only on those terms. So,
U = (β^t)log((Subscript[w, t])+(1+(Subscript[r, t])-δ)(Subscript[K, t])+(Subscript[π, t])-(Subscript[K, t+1])) + (β^(t+1))log((Subscript[w, t+1])+(1+(Subscript[r, t+1])-δ)(Subscript[K, t+1])+(Subscript[π, t+1])-(Subscript[K, t+2]))
Now we can differentiate with respect to (Subscript[K, t+1]) and set equal to 0.

U = (β^t)*Log[(Subscript[w, t])+(1+(Subscript[r, t])-δ)*(Subscript[K, t])+(Subscript[π, t])-(Subscript[K, t+1])]+(β^(t+1))*Log[(Subscript[w, t+1])+(1+(Subscript[r, t+1])-δ)*(Subscript[K, t+1])+(Subscript[π, t+1])-(Subscript[K, t+2])];
D[U,(Subscript[K, t+1])]==0

The output from differentiation above is:
					
We can rearrange to get:
((β^t)/((Subscript[w, t])+(1+(Subscript[r, t])-δ)(Subscript[K, t])+(Subscript[π, t])-(Subscript[K, t+1]))) = (((β^(t+1))(1+(Subscript[r, t+1])-δ))/((Subscript[w, t+1])+(1+(Subscript[r, t+1])-δ)*(Subscript[K, t+1])+(Subscript[π, t+1])-(Subscript[K, t+2])))
Then we divide both sides by (β^t) to get:
(1/((Subscript[w, t])+(1+(Subscript[r, t])-δ)(Subscript[K, t])+(Subscript[π, t])-(Subscript[K, t+1]))) = ((β(1+(Subscript[r, t+1])-δ))/((Subscript[w, t+1])+(1+(Subscript[r, t+1])-δ)*(Subscript[K, t+1])+(Subscript[π, t+1])-(Subscript[K, t+2])))
Then we can realise that the denominator of the LHS is (Subscript[C, t]) and the denominator of the RHS is (Subscript[C, t+1]) so,
(1/(Subscript[C, t])) = ((β(1+(Subscript[r, t+1])-δ))/( (Subscript[C, t+1])))
Then multiply both sides by (Subscript[C, t+1]) to get:
((Subscript[C, t+1])/(Subscript[C, t])) = β(1+(Subscript[r, t+1])-δ)
This gives us the Euler equation

To find the FOC for the firm we take the profit maximising condition and substitute (Subscript[Y, t]) into it. So we get:
(Subscript[π, t]^f) = (max_((Subscript[L, t]^f),(Subscript[K, t]^f)))  ((((Subscript[K, t]^f)))^α)((((Subscript[AL, t]^f)))^(1-α))- (Subscript[w, t])(Subscript[L, t]^f) - (Subscript[r, t])(Subscript[K, t]^f)
Then we can differentiate with respect to (Subscript[L, t]^f) and (Subscript[K, t]^f) to get the first order conditions.

Profit=((((Subscript[K, t])))^α)(((A*(Subscript[L, t])))^(1-α))-(Subscript[w, t])(Subscript[L, t])-(Subscript[r, t])(Subscript[K, t]);
D[Profit,(Subscript[L, t])]==0
D[Profit,(Subscript[K, t])]==0

The output from differentiation in Mathematica gives us:
						
We can rearrange this to get:
(Subscript[w, t]) = A (1-α) ((((Subscript[K, t]^f)))^α) ((((Subscript[AL, t]^f)))^(-α)) = A (1-α) ((((Subscript[K, t]^f)))^α) ((((1/(Subscript[AL, t]^f))))^α) = A (1-α) (((((Subscript[K, t]^f)/(Subscript[AL, t]^f))))^α)
(Subscript[r, t]) = α ((((Subscript[K, t]^f)))^(-1+α)) ((((Subscript[AL, t]^f)))^(1-α)) = α ((((1/(Subscript[K, t]^f))))^(1-α)) ((((Subscript[AL, t]^f)))^(1-α)) = α (((((Subscript[AL, t]^f)/(Subscript[K, t]^f))))^(1-α))

Part C. The solution to the growth model we studied in class was characterized
by the Euler equation written only in terms of capital stocks. Write out this Euler
equation.

From above we know that the Euler equation is : ((Subscript[C, t+1])/(Subscript[C, t])) = β(1+(Subscript[r, t+1])-δ).
We know that (Subscript[Y, t]) = (Subscript[C, t]) + (Subscript[I, t]),  (Subscript[I, t])= (Subscript[K, t+1])-(1-δ)(Subscript[K, t]),  (Subscript[r, t]) = α (((((Subscript[AL, t])/(Subscript[K, t]))))^(1-α)) and  (Subscript[r, t+1]) = α (((((Subscript[AL, t+1])/(Subscript[K, t+1]))))^(1-α)).
(Subscript[C, t]) = (Subscript[Y, t]) - (Subscript[I, t]) = ((((Subscript[K, t])))^α)((((Subscript[AL, t])))^(1-α)) - ((Subscript[K, t+1])-(1-δ)(Subscript[K, t]))
Similarly, for (Subscript[C, t+1]):
(Subscript[C, t+1]) = ((((Subscript[K, t+1])))^α)((((Subscript[AL, t+1])))^(1-α)) - ((Subscript[K, t+2])-(1-δ)(Subscript[K, t+1]))
Substituting this into the Euler equation:
((((((Subscript[K, t+1])))^α)((((Subscript[AL, t+1])))^(1-α))-((Subscript[K, t+2])-(1-δ)(Subscript[K, t+1])))/(((((Subscript[K, t])))^α)((((Subscript[AL, t])))^(1-α))-((Subscript[K, t+1])-(1-δ)(Subscript[K, t])))) = β(1+α((((A/(Subscript[K, t+1]))))^(1-α))-δ)
We know that L is equal to 1 so:
((((((Subscript[K, t+1])))^α)(((A))^(1-α))-(Subscript[K, t+2])+(1-δ)(Subscript[K, t+1]))/(((((Subscript[K, t])))^α)(((A))^(1-α))-(Subscript[K, t+1])+(1-δ)(Subscript[K, t]))) = β(1+α((((A/(Subscript[K, t+1]))))^(1-α))-δ)
which is the Euler equation written only in terms of capital stock.

Part d. Recall, a steady state of the above problem was an equilibrium where all variables were constant. Use the first order condition above to solve for the steady state level of capital (SuperStar[K])

A steady state is a competitive equilibrium where all variables are constant:
(Subscript[K, t+1]) = (Subscript[K, t]) = (SuperStar[K])
(Subscript[C, t+1]) = (Subscript[C, t]) = (SuperStar[C])
(Subscript[r, t+1]) = (Subscript[r, t]) = (SuperStar[r])
(Subscript[w, t+1]) = (Subscript[w, t]) = (SuperStar[w])
From the Euler equation in steady state it is true that:
((Subscript[C, t+1])/(Subscript[C, t])) = β(1+(Subscript[r, t+1])-δ) --> ((SuperStar[C])/(SuperStar[C])) = β(1+(SuperStar[r])-δ) hence,
1 = β(1 + (r^( * )) - δ) Divide both sides by β
(1/β) = 1 + (SuperStar[r]) - δ  Move  (SuperStar[r]) to the LHS and (1/β) to the RHS
(SuperStar[r]) = (1/β) - 1 + δ
From the first order condition we know:
(Subscript[r, t]) = α ((((A/(Subscript[K, t]))))^(1-α)) --> (SuperStar[r]) = α ((((A/(SuperStar[K]))))^(1-α))
This implies that in the steady state:
(1/β) - 1 + δ = α ((((A/(SuperStar[K]))))^(1-α)) Divide both sides by α
(((1/β)-1+δ)/α) = ((((A/(SuperStar[K]))))^(1-α))  Expand the RHS
(((1/β)-1+δ)/α) = ((A^(1-α))/((((SuperStar[K])))^(1-α)))  Multiply both sides by ((((SuperStar[K])))^(1-α))
((((SuperStar[K])))^(1-α))(((1/β)-1+δ)/α) = (A^(1-α))  Multiply both sides by (α/((1/β)-1+δ))
((((SuperStar[K])))^(1-α)) = (A^(1-α)) ((α/((1/β)-1+δ)))  Take both sides to the power of (1/(1-α))
(SuperStar[K]) = A  ((((α/((1/β)-1+δ))))^(1/(1-α)))  We have found the steady state level of capital (SuperStar[K])

Part e. Assume that the model converges to the steady state after T=70
periods.Write code that will allow you to solve the above model using the shooting
algorithm that we studied in class. For the parameters assume α=0.33,β=0.98,δ=0.1,A=1 and K0=0.1.Plot the evolution of the capital stock you find.

α=0.33;β=0.98;δ=0.1;A=1;
H[t_]:=(((1-δ) (Subscript[k, 1+t])+(Subscript[k, 1+t]^α)(A^(1-α))-(Subscript[k, 2+t]))/((1-δ) (Subscript[k, t])+(Subscript[k, t]^α)(A^(1-α))-(Subscript[k, 1+t])))-β (1-δ+α ((((A/(Subscript[k, 1+t]))))^(1-α)) );
T=70;
(Subscript[k, 0])=0.1;(Subscript[k, T+2])=(Subscript[k, T+1]);
InitialGuess=Table[{(Subscript[k, i]),0.2},{i,1,T+1}];
answerk=FindRoot[Table[H[t]==0,{t,0,T}],InitialGuess];
answerkvector=Table[(Subscript[k, i]),{i,1,Length[answerk]}]/.answerk;
ListPlot[{answerkvector},Joined->True,PlotRange->All]
answerkvector

Genetic Algorithm

Part a. Suppose that you wish to solve the above model using a genetic algorithm. Assume that the model converges to the steady state after T periods.
Create a module called score that takes a potential solution (as well as T) and
returns the fitness of that particular solution. Do not forget that this function
must account for the fact that you want your model to converge after T periods!

(*Define constants*)
A=1; (*Technology parameter*)
δ=0.1; (*Depreciation rate*)
α=0.33; (*Capital share in production function*)
K0=0.1; (*Initial capital stock*)
β=0.98; (*Discount factor*)

(*Module to compute capital,output,consumption,and investment over time*)
answers[S_List,T_]:=Module[{K,Y,s,I,C,YList,CList,KList,IList,i},(*Initialize lists to store results over time*)
KList={K0};(*Start with initial capital*)
YList={};(*List for output over time*)
CList={};(*List for consumption over time*)
IList={};(*List for investment over time*)
K=K0;(*Initialize capital*)
(*Simulation loop for T periods*)
For[i=0,i<=T,i++,
Y=K^α A^(1-α);(*Compute output (Y_t) using the production function*);
AppendTo[YList,Y];(*Store output in the list*)
I=S[[i+1]] Y;(*Compute investment (I_t) as a share of output*)
AppendTo[IList,I];(*Store investment in the list*);C=(1-S[[i+1]]) Y;(*Compute consumption (C_t) as the remaining share of output*);
AppendTo[CList,C];(*Store consumption in the list*)
K=(1-δ) K+I;(*Update capital for the next period*)
AppendTo[KList,K]; (*Store updated capital in the list*)];
{KList,YList,CList,IList}(*Return all computed lists*)];

(*Module to compute the total fitness score*)
Score[S_List,T_]:=Module[{KList,YList,CList,IList,i,fitness,utility,Totalfitness},(*Call the answers module to get computed lists*){KList,YList,CList,IList}=answers[S,T];
fitness={};(*Initialize the fitness list to store utilities over time*);
(*Loop to calculate discounted utility for each period*)
For[i=0,i<=T,i++,
utility=β^i Log[CList[[i+1]]];(*Compute discounted utility*)
AppendTo[fitness,utility]; (*Store utility in the fitness list*)];
(*Compute total fitness as the sum of utilities+terminal utility*)Totalfitness=Total[fitness]+β^T ((Log[CList[[T+1]]])/(1-β));
Return[Totalfitness];(*Return the total fitness score*)];

Part b .  Assume that the model converges to the steady state after T = 70
periods. Write code that will allow you to solve the above model using the genetic
algorithm that we studied in class. Be sure to convert your potential solutions
to binary and apply mutations on the binary code. Your code should contain (at
least) the following modules (but may contain others):
• selectionGen - chooses which chromosomes are selected and which are dropped
• parentsGen - chooses which selected chromosomes become parents, given their
fitness scores
• crossover - creates child chromosomes from selected parent chromosomes
• mutateChromosome - mutates the child chromosomes
Run your code for at least a 1000 generations. Choose your population size to be
at least 70. Choose a mutation rate of 2.5%. For the parameters assume α = 0.33,
β = 0.98, δ = 0.1, A = 1 and K0 = 0.1. Once solved, plot the evolution of the
capital stock. ALSO (on another graph) plot the evolution of the fitness function.

(*Function to perform selection based on fitness scores*)selectionGen[pop0_,fitness0_]:=Module[{pop=pop0,fitness=fitness0,location,survpop},(*Determine the indices of the top 50% of the population by fitness*)location=Ordering[fitness,-(Length[pop]/2)];
(*Select the individuals at those indices to form the surviving population*)survpop=pop[[location]];
Return[survpop];(*Return the surviving population*)];

(*Function to generate parent pairs for reproduction based on fitness*)parentsGen[survpop0_]:=Module[{survpop=survpop0,Probabilty,parentpop },(*Compute selection probabilities for each individual in the surviving population*)Probabilty=1.0*Table[Length[survpop]-i+1,{i,1,Length[survpop]}]/Sum[i,{i,1,Length[survpop]}];(*Select parent pairs using the computed probabilities*)(*RandomChoice selects individuals with weighted probabilities*)(*{Length[survpop]/2,2} creates pairs of parents for half the population*)parentpop=RandomChoice[Probabilty->survpop,{Length[survpop]/2,2} (*Assign probabilities to each individual and generate pairs of parents*)];
Return[parentpop];(*Return the selected parent pairs*)];

(*Function to perform crossover on a list of parent pairs*)crossover[parents_]:=Module[{numberToBinary,crossoverBinary,results},(*Function to convert a number to binary with a specified number of bits*)numberToBinary[num_]:=IntegerDigits[Floor[num*2^21],2,21];
(*Function to perform crossover between two binary parents*)crossoverBinary[parent1_,parent2_]:=Module[{crossoverPoint,child1,child2},(*Randomly choose a crossover point within the binary sequence*)crossoverPoint=RandomInteger[{0,Length[parent1]}];
(*Create two offspring by swapping parts of the binary parents*)child1=Join[parent1[[1;;crossoverPoint]],parent2[[crossoverPoint+1;;]]];
child2=Join[parent2[[1;;crossoverPoint]],parent1[[crossoverPoint+1;;]]];
{child1,child2}(*Return the two offspring*)];
(*Process each pair of parents in the list*)results=Table[Module[{binaryList1,binaryList2  },(*Convert numerical values of the current pair of parents to binary*)binaryList1=numberToBinary/@parents[[i,1]];
binaryList2=numberToBinary/@parents[[i,2]];
(*Perform crossover for all corresponding binary sequences*)Table[crossoverBinary[binaryList1[[j]],binaryList2[[j]]],{j,Length[binaryList1]}]],{i,Length[parents]} (*Iterate over all pairs of parents*)];
Return[results];(*Return the results of the crossover process*)];

(*Function to perform mutation on a list of chromosomes*)mutateChromosome[chromosomes_,mutationRate_]:=Module[{mutateBit,mutateList,MC,binaryToDecimal },(*Function to randomly mutate a single bit*)mutateBit[bit_]:=If[RandomReal[]<mutationRate,1-bit,bit (*Compare a random value with the mutation rate*)(*Flip the bit if mutation occurs (0 becomes 1,1 becomes 0)*)(*Keep the bit unchanged otherwise*)];
(*Function to apply mutation recursively to all binary lists*)mutateList[list_]:=If[Head[list]===List,(*Check if the current element is a list*)mutateList/@list,(*If so,apply mutation recursively to each element*)mutateBit[list] (*Otherwise,apply mutation to the bit*)];
(*Apply mutation to all chromosomes*)
MC=Flatten[mutateList[chromosomes],2];(*Flatten the mutated lists to a depth of 2*)(*Function to convert a binary list to a decimal number*)binaryToDecimal[binary_List]:=Total[Table[binary[[i]]*2^(-i),(*Convert binary to decimal using positional weights*){i,1,Length[binary]} (*Iterate over all bits in the binary list*)]];
(*Convert the mutated binary chromosomes back to decimal form*)Partition[binaryToDecimal/@MC,71]];

mutationRate=0.025; (*Probability of mutation for each bit*)
T=70; (*Number of time periods untill steady state*)
Population=Table[RandomReal[{0.00001,0.99999},T+1],80];(*Initialize population*)

(*Initialize memory to track progress*)
fitnessMemoryMean={}; (*To store the mean fitness over generations*)
fitnessMemoryMax={}; (*To store the maximum fitness over generations*)

(*Evolutionary algorithm:Loop over generations*)
For[gen=1,gen<=100,gen++,(*Run for 1500 generations*)
fitness=Score[#,T]&/@Population;(*Calculate fitness for the current population*);(*Compute fitness for each individual*);
fittest=Population[[Ordering[fitness,1][[1]]]];(*Identify the fittest individual in the population*)(*Track mean and max fitness values over generations*)fitnessMemoryMean=Append[fitnessMemoryMean,Mean[fitness]];fitnessMemoryMax=Append[fitnessMemoryMax,Max[fitness]];
popSurvivors=selectionGen[Population,fitness];(*Select survivors based on fitness*);popParents=parentsGen[popSurvivors];(*Generate parents for crossover*);preMutationChildren=crossover[popParents];(*Perform crossover to generate children*);(*Apply mutation to the children*)postMutationChildren=1.0*mutateChromosome[preMutationChildren,mutationRate];
(*Update population by combining survivors and mutated children*)Population=Join[popSurvivors,postMutationChildren];];

(*Final fitness evaluation after all generations*)
fitness=Score[#,T]&/@Population; (*Calculate final fitness for the population*)

(*Track final mean and maximum fitness*)
fitnessMemoryMean=Append[fitnessMemoryMean,Mean[fitness]]; (*Final mean fitness*)
fitnessMemoryMax=Append[fitnessMemoryMax,Max[fitness]]; (*Final max fitness*)

(*Retrieve the fittest individual in the final population*)
fittestanswer=Population[[Ordering[fitness,1][[1]]]]; (*Best individual*)

(*Return the fittest individual as the result*)
fittestanswer

{KList,YList,CList,IList}=answers[fittestanswer,T];
KList;(*Capital stock at last generation*)
ListPlot[{Legended[KList,"Genetic Capital Stock"]},Joined->True,PlotRange->All](*Evolution of capital stock*)
ListPlot[{Legended[fitnessMemoryMax,"Maximum Fitness"],Legended[fitnessMemoryMean,"Mean Fitness"]},Joined->True,PlotRange->All](*Evolution of the fitness function*)

Part C.  Plot the evolution of the capital stock from the shooting - model and the
capital stock from the genetic algorithm over time on the same graph . Also include
a line to indicate the capital stock in the steady state . On another graph plot the
evolution of both the savings rate from the shooting - model and the savings rate
from the genetic algorithm . What do you notice?

ListPlot[{Legended[KList,"Genetic Capital Stock"],Legended[Prepend[answerkvector,K0],"Shooting Capital Stock"]},Joined->True,Epilog->{Dashed,Black,Line[{{0,A  ((((α/((1/β)-1+δ))))^(1/(1-α)))},{71,A  ((((α/((1/β)-1+δ))))^(1/(1-α)))}}]}]
(*Evolution of capital stock from shooting algorithm and genetic algorithm*)

(*Calculations to get savings rate for the shooting algorithm using our answerkvector*)
yList=((A^(1-α))*Prepend[answerkvector,K0]^α);
IList={};
For[t=1,t<72,t++,
AppendTo[IList,Prepend[answerkvector,K0][[t+1]]-(1-δ)* Prepend[answerkvector,K0][[t]]];];
yList1=Drop[yList,{71}];
SList=IList/yList1;
ListPlot[{Legended[fittestanswer,"Genetic Savings Rate"],Legended[SList,"Shooting Savings Rate"]},Joined->True,PlotRange->{0,1}](*Evolution of savings rate from shooting model and genetic algorithm*)

Capital Stock Evolution : The capital stock from the shooting model (the smooth orange curve) converges steadily toward the steady state . The capital stock from the genetic algorithm (the blue curve) fluctuates more significantly and does not converge as smoothly . These fluctuations likely arise from the stochastic nature of the genetic algorithm, which introduces randomness via crossover and mutation . 
Savings Rate Evolution : The savings rate from the shooting model (orange curve) declines smoothly over time, reflecting optimal behavior determined analytically . The savings rate from the genetic algorithm (blue curve) shows much greater volatility and does not follow the smooth path of the shooting model . This again highlights the trade - off in genetic algorithms between optimization and exploration, with the randomness causing deviations from the theoretical path . Insights : The shooting model provides a smooth and predictable trajectory for both capital stock and savings rate, which aligns with theoretical expectations . The genetic algorithm introduces variability due to its probabilistic operations . While it eventually converges toward reasonable values, the fluctuations in both capital stock and savings rate highlight the less precise nature of this method . Including the steady - state line provides a clear visual reference for assessing the convergence of both approaches . The genetic algorithm, while approximating the steady state, does so with some degree of variability, unlike the shooting model, which reaches the steady state more precisely . 
Suggestions for Improvement : To reduce the fluctuations in the genetic algorithm' s results, consider : Increasing the population size . Running the algorithm for more generations . Adjusting the mutation rate (e . g ., lowering it slightly) .