れん９８７６ - 経済統計練習帳

beta = 0.95; 
a = 0.33; 
z1 = 1/ beta - 1; 
ks = 0.5*pow(z1/a,1/(a-1)); 
h = 2 * ks / 100; 
for ( n=1; n<=100; n++ ){ 
k[n] = n * h; 
lx[n] = 0.5; 
cx[n] = pow(k[n] , a) * pow(lx[n] , 1 - a); 
}   
t1 = 0; 
do { 
for ( n=10; n<=90; n++ ){ 
k1 = k[n] + pow(k[n] , a) * pow(lx[n] ,1 - a) - cx[n]; 
n1 = k1 / h; 
n2 = floor(n1); 
n3 = n2 + 1; 
c1 = cx[n2] + (n1 - n2) * (cx[n3] - cx[n2]); 
l1 = lx[n2] + (n1 - n2) * (lx[n3] - lx[n2]); 
r1 = a * pow(k1 , a - 1) * pow(l1,1 - a); 
cp[n] = c1 / (beta * (1 + r1)) ; 
w1 = (1 - a) * pow(k[n] , a) * pow(lx[n] , -a); 
lp[n] = 1 - cx[n] / w1; 
} 
e = 0; 
for ( n=10; n<=90; n++ ){ 
e = e + pow(cx[n] - cp[n],2) + pow(lx[n] - lp[n],2); 
} 
for ( n=10; n<=90; n++ ){ 
cx[n] = cp[n]; 
lx[n] = lp[n]; 
} 
If ( e < 0.0001) {t1 = 1000; 
} 
t1 = t1 + 1; 
} 
while(t1 < 100);   
ms = 10; 
for ( n=10; n<=90; n++ ){ 
px[n] = 1; 
} 
t2 = 0; 
do { 
for ( n=10; n<=90; n++ ){ 
p1 = 0.9 * px[n]; 
p2 = 1.1 * px[n]; 
k1 = k[n] + pow(k[n] , a) * pow(lx[n] ,1 - a) - cx[n]; 
n1 = k1 / h; 
n2 = floor(n1); 
n3 = n2 + 1; 
pz = px[n2] + (n1 - n2) * (px[n3] - px[n2]); 
pi = pz / p1 - 1; 
l1 = lx[n2] + (n1 - n2) * (lx[n3] - lx[n2]); 
r1 = a * pow(k1 , a - 1) * pow(l1,1 - a); 
ix = (1 + r1) * (1 + pi) - 1; 
z1 = ms * ix / (cx[n] * (1 + ix)) - p1; 
t3 = 0; 
do{ 
pi = pz / p2 - 1; 
ix = (1 + r1) * (1 + pi) - 1; 
z2 = ms * ix / (cx[n] * (1 + ix)) - p2; 
p3 = p2 - z2 * (p2 - p1) / (z2 - z1); 
p1 = p2; 
p2 = p3; 
z1 = z2; 
If (pow(z2 , 2) < 0.0001) { 
t3 = 1000; 
} 
t3 = t3 + 1; 
} 
while(t3 < 100); 
ps[n] = p2; 
} 
e = 0; 
for ( n=10; n<=90; n++ ){ 
e = e + pow(px[n] - ps[n], 2); 
} 
for ( n=10; n<=90; n++ ){ 
px[n] = ps[n]; 
} 
If ( e < 0.0001){ 
t2 = 1000; 
} 
t2 = t2 + 1; 
} 
while(t2 < 100); 
t=0; 
kt[t]=k[30]; 
for ( t=0; t<=99; t++ ){ 
n1 = kt[t] / h; 
n2 = floor(n1); 
n3 = n2 + 1 ; 
ct[t] = cx[n2] + (n1 - n2) * (cx[n3] - cx[n2]); 
lt[t] = lx[n2] + (n1 - n2) * (lx[n3] - lx[n2]); 
pt[t] = px[n2] + (n1 - n2) * (px[n3] - px[n2]); 
kt[t+1] = kt[t] + pow(kt[t] , a) * pow(lt[t] ,1 - a) - ct[t]; 
} 
for ( t=1; t<=98; t++ ){ 
print(pt[t]); 
print(","); 
} 
t=99; 
print(pt[t]);