#Run iterations to solve for the initial equilibrium values

for (j in 1:50){

for (t in 2:(T-1)){

K[t]<-s[t-1]*L[t-1]
k[t]<-K[t]/L[t]

Y[t]<-K[t]^alpha*L[t]^(1-alpha)
r[t]<-alpha*k[t]^(alpha-1)-1
w[t]<-(1-alpha)*k[t]^alpha

s[t]<-(w[t]*(1-tau[t])-(1+rho)/(1+r[t+1])*w[t+1]*tau[t+1]*(1+n[t+1]))/(2+rho)

C_y[t]<-(1-tau[t])*w[t]-s[t]
C_o[t]<-(1+r[t])*s[t-1]+tau[t]*w[t]*(1+n[t])

U[t]<-log(C_y[t])+log(C_o[t+1])/(1+rho)

}

r[T]<-r[T-1]
w[T]<-w[T-1]

s[1]<-s[T-1]

plot(U[95:110], col=j, type="l")
}