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

for ( m=1; m<=10; m++ ){ 
th[m] = 0.95 + 0.01 * m; 
} 
beta = 0.95; 
a = 0.33; 
ks = pow((1 / beta - 1) / a , 1 / (a - 1)); 
h = 2 * ks / 100; 
for ( n=1; n<=100; n++ ){ 
k[n] = n * h; 
} 
for ( m=1; m<=10; m++ ){ 
for ( n=1; n<=100; n++ ){ 
cx[m][n] = th[m] * pow(k[n] ,a); 
} 
} 
t1 = 0; 
Do { 
for ( m=1; m<=10; m++ ){ 
for ( n=10; n<=90; n++ ){ 
uc = 0; 
for ( s=1; s<=10; s++ ){ 
k1 = k[n] + th[m] *pow(k[n] ,a) - cx[m][n]; 
r1 = th[s] * a * pow(k1 , a - 1) ; 
n1 = k1 / h; 
n2 = floor(n1); 
n3 = n2 + 1; 
c1 = cx[s][n2] + (n1 - n2) * (cx[s][n3] - cx[s][n2]); 
uc = uc + (beta * (1 + r1)) / c1; 
} 
uc = uc / 10; 
cp[m][n] = 1 / uc; 
} 
} 
e = 0; 
for ( m=1; m<=10; m++ ){ 
for ( n=10; n<=90; n++ ){ 
e = e + pow(cx[m][n] - cp[m][n], 2); 
} 
} 
for ( m=1; m<=10; m++ ){ 
for ( n=10; n<=90; n++ ){ 
cx[m][n] = cp[m][n]; 
} 
} 
If (e < 0.0001) { 
t1 = 1000; 
} 
t1 = t1 + 1; 
} 
while(t1 < 100); 
 
ms = 50; 
for ( m=1; m<=10; m++ ){ 
for ( n=10; n<=90; n++ ){ 
px[m][n] = 1; 
} 
} 
 
t2 = 0; 
Do { 
for ( m=1; m<=10; m++ ){ 
for ( n=10; n<=90; n++ ){ 
p1 = 1.01 * px[m][n]; 
p2 = 0.99 * px[m][n]; 
k1 = k[n] + th[m] *pow(k[n] ,a) - cx[m][n]; 
n1 = k1 / h; 
n2 = floor(n1); 
n3 = n2 + 1; 
 
um = 0; 
for ( s=1; s<=10; s++ ){ 
r1 = a * th[s] * pow(k1 ,a - 1); 
c1 = cx[s] [n2] + (n1 - n2) * (cx[s][n3] - cx[s][n2]); 
pc = px[s] [n2]+ (n1 - n2) * (px[s][n3] - px[s][n2]); 
pp = pc / p1; 
i = pp * (1 + r1) - 1 ; 
um = um + c1 / (beta * i * pp); 
} 
um = 0.1 * um; 
 
z1 = ms / um - p1; 
 
t3 = 0; 
Do { 
 
um = 0; 
for ( s=1; s<=10; s++ ){ 
r1 = a * th[s] * pow(k1 , a - 1); 
c1 = cx[s] [n2] + (n1 - n2) * (cx[s][n3] - cx[s][n2]); 
pc = px[s] [n2]+ (n1 - n2) * (px[s][n3] - px[s][n2]); 
pp = pc / p2; 
i = pp * (1 + r1) - 1;   
um = um + c1 / (beta * i * pp); 
} 
um = 0.1 * um; 
 
z2 = ms / um - p2; 
p3 = p2 - z2 * (p2 - p1) / (z2 - z1); 
p1 = p2; 
p2 = p3; 
z1 = z2; 
If (pow(z1, 2) < 0.01){ 
t3 = 1000; 
} 
t3 = t3 + 1; 
} 
while (t3 < 100); 
ps[m][n] = p2; 
} 
} 
 
e = 0; 
for ( m=1; m<=10; m++ ){ 
for ( n=10; n<=90; n++ ){ 
e = e + pow(px[m][n] - ps[m][n],2); 
} 
} 
for ( m=1; m<=10; m++ ){ 
for ( n=10; n<=90; n++ ){ 
px[m][n] = ps[m][n]; 
} 
} 
If (e < 0.01){ 
t2 = 1000; 
}   
t2 = t2 + 1; 
} 
while ( t2 < 100); 
t=1; 
kt[t] = k[50]; 
for ( t=1; t<=100; t++ ){ 
n1 = kt[t] / h; 
n2 = floor(n1); 
n3 = n2 + 1 ; 
m = rand(1,10); 
tht[t] = th[m]; 
ct[t] = cx[m][n2] + (n1 - n2) * (cx[m][n3] - cx[m][n2]); 
pt[t] = px[m][n2] + (n1 - n2) * (px[m][n3] - px[m][n2]); 
kt[t+1] = kt[t] + tht[t]*pow(kt[t] , a) - ct[t]; 
}