ぃね - 経済統計練習帳

import numpy as np
import math
from matplotlib import pyplot
a=0.33
beta=0.95
ks=(1/beta - 1)/a
q=1/(a-1)
ks=ks**q
h=ks/100
k=np.zeros(101)
for i in range(101):
    k[i]=i*h
cx=np.zeros(101) 
cxs=np.zeros(101)
for i in range(101):
        cx[i]=k[i]**a
 
for tr in range(25):
 
    for i in range(10,91):
        c1=cx[i]
        k1=k[i]+k[i]**a-c1
        n1=k1/h
        if n1>90:n1=90
        n2=math.floor(n1)
        n3=n1-n2        
        c2=cx[n2]+n3*(cx[n2+1]-cx[n2])
        r1=a*k1**(a-1)
        uc=(beta*(1+r1))/c2
        cxs[i]=1/uc
    e=0
    for i in range(101):
        e=e+(cx[i]-cxs[i])*(cx[i]-cxs[i])
    print(e)
 
    for i in range(101):
        cx[i]=cxs[i]
 
 
pyplot.plot(k, cx)
pyplot.show()