Numeric Recipes for Econometrics(0) SHIMURA Masato jcd02773@nifty.com 2010 12 2 1 3-2 1.1...................................... 2 1.2 3........................... 8 1.3....................................... 14 1.4 AR.................................... 16 1.5.......................................... 19 2 20 2.1........................................ 20 2.2...................................... 21 2.3..................................... 25 2.4...................................... 27 2.5...................................... 29 2.6.......................................... 30 2.7.................................. 31 2.8................................... 34 2.9........................................ 35 2.10.......................................... 36 A 38 B EXCEL 39 B.1............................................ 39 1
B.2 EXCEL............................... 39 B.3 tara............................................. 39 B.4 CSV....................................... 41. 1 3 - (AR) y X X X 1.1 1.1.1 X 1 X 2 X 3 X n,.y (y Example (i js) DN11 *1 5.6 30 Stich(,.) 5.8 26 (,) 6 33 EX0=: 1 2 3,4 5 6,:7 8 9 6.2 31 6.4 33 ijs =: 6.4 35 6.4 37 6.6 36 6.8 33 EXCEL CSV *2 EXCEL CSV APPENDIX ( ) (GDP) () () () 94 473 272 65 131 95 484 278 67 138 *1 = (equal) *2 Camma Separared Value 2
96 502 286 74 142 97 504 283 79 137 98 500 285 75 134 99 503 285 75 137 0 515 287 81 141 1 509 289 78 137 2 514 291 76 139 3 524 293 82 137 4 534 296 86 137 *3 ( ) 2 2006 CSV *4 require files csv plot ] DN10=. ".@> readcsv /data/excel/stat_j/csv/ban_1.csv /data/excel/stat_j/csv/ 1.1.2 EXCEL CSV =. (L:0) DN10 = ;( 1) ( ;( 2)) = ".@>. J 1 2 (J J, 0 ( 1) *3 numeric recipe data.ijs *4 CSV Camma Separated Value () EXCEL DB CSV 3
( Take { Take f irst {. 0{ 1 DN11 Tale last {: {: 1 DN11 Drop f irst }. }. 1 DN11 Drop last }: } : 1 DN11 J n * 5 f unction name U sage 0 puc pickup column rec remove column puc set pickup with set 1 puc1 pickup column rec1 remove column puc set1 pickup with set puc=:pick_column=: [ {"1 ] puc1=:pick_column=: (>:@[) {"1 ] rec=:remove_column=: 4 : (I.-.(i.{:@$ y ) e. x ) puc y NB. remove_raw rer=: remove_raw=: 4 : (I.-.(i. # y ) e. x ) { y Example ] a=.?. 7 4 $ 28 6 3 7 8 2 23 0 25 20 26 18 14 13 6 3 0 2 4 26 6 19 16 8 16 8 11 5 27 2 3 puc a 7 8 0 25 18 14 3 0 26 6 8 16 5 27 1 2 rec a 6 8 2 25 20 14 13 0 2 6 19 16 8 27 *5 0, 1, 2, 3 0 1, 2, 3.. 1 0 1 4
y X y 2 2 1,.X 0,.X 1,.X 2,.X n ( 1 ) y %. 1,.X 1,.X,.y y reg_ols=: %. 1&,.@] y, X J x y type 1 type 2 y is 0, x is 1 (0 puc a) reg_ols 1 puc a 5.72402 0.336313 f = 5.72402 + 0.336313x y is 0, x is 1 2 3 ( ) (0 puc a) reg_ols 0 rec a 8.92071 0.597347 _0.065022 _0.429712 y = 8.92071 + 0.597347x 0 0.065022x 1 0.429712x 2 y is 0 x is 3 (0;3) puc_set a x y 3 0 puc a 8 6 x y 25 2 8 6 14 20 25 2 0 13 14 20 6 2 0 13 16 19 6 2 27 8 16 19 regx (0;3) puc_set a 27 8 11.1634 _0.0848279 y = 11.1634 0.0848279x y is 0 x is 1 3 (0; 1 3) puc_set a x0 x1 y 5
3 8 6 23 25 2 26 14 20 6 0 13 4 6 2 16 16 19 11 27 8 regx (0;1 3) puc_set a 8.05491 0.588484 _0.403745 f = 8.05491 + 0.588484x 0 0.403745x 1 1.1.3 *6 EPS EMF J sin 1 0.8 0.6 0.4 0.2 0-0.2-0.4-0.6-0.8-1 -5-4 -3-2 -1 0 1 2 3 4 5 title sin plot _5 5 ; 1&o. NB. sin from _5 to 5 pd eps /temp/sin0.eps NB. save by eps 5 5 plot (/temp ) J Grammar *6 J 6
". format ": ".@> >open Box < Box require readcsv csv require files csv plot pd plot driver eps 7
1.2 3 1.2.1 Matrix divide (AR) 3 %. matrix Divide K.E.Iverson *7, Ordinaly Least regx=: 3 : ({:"1 y ) %. 1,.}:"1 y Square OLS Polynomial poly1=:4 : y %. (>:i. # y )ˆ/i. >: x Auto Regression AR ar0=:4 : (x }.tmp) %..("1)}: >x <\ tmp=:y -(+/%#)y 1.2.2 OLS y 1 1 x 11 y 1 1 x 11 x 21 x k1 y = y 2 y 3 y n, 1 x 12 1 x 13.. 1 x 1n y = y 2 y 3 y n, 1 x 12 x 22 x k2 1 x 13 x 23 x k3....... 1 x 1n x 2n x kn (X X)ˆβ = X y *7 APL /. 8
ˆβ = (X X) 1 X y (X X) 1 X y = X y X X = y X X 1 Working Example DN11 X cm Y (%) y 1,.X DN11 K.E.Iverson matrix divide (%.) 9
1.2.3 X,.Y X 1 y X Y 1,.X %. DN11 X Y ------- 5.6 30 5.8 26 6 33 6.2 31 6.4 33 6.4 35 6.4 37 6.6 36 6.8 33 1,.0 puc DN11 1 5.6 1 5.8 1 6 1 6.2 1 6.4 1 6.4 1 6.4 1 6.6 1 6.8 y %. 1,.X (1 puc DN11) %. 1,. 0 puc DN11 _5.01128 6.03383 reg0 DN11 _5.01128 6.03383 y = 5.01128 + 6.03383x *8 Script reg0 reg0=:3 : 0 NB. select trend d ata or multi data if. 1= +/ * $ y do. reg_t y elseif. do. regx y end. ) reg_t=:3 : y %. 1,. >: i. # y regx=:3 : ({:"1 y ) %. 1,.}:"1 y 1.2.4 *8 X Y reg_t reg0 DN11 regx 10
X 1 (X X) 1 X y = X y X X = y X %. ({:"1 DN11) %. {."1 DN11 5.23401 (1 puc DN11) %. 0 puc DN11 5.23401 y = 5.23401x 36 34 32 30 28 26 5.6 5.8 6 6.2 6.4 6.6 6.8 1 1.2.5 reg exam ad reg0 reg_exam_ad *9 reg0 reg_exam_ad DN11 +------+----------------+ f= _5.01128 6.03383 *9 ( ) J 11
+------+----------------+ corr=: 47.8237 +------+----------------+ AIC: 18.8684 AIC +------+----------------+ DW= 2.17534 / +------+----------------+ t=: _0.33633 2.53299 t +------+----------------+ 1.2.6 line f it reg0 linefit_reg0 DN11 pd eps /temp/reg_0.eps NB. Save 36 34 32 30 28 26 5.6 5.8 6 6.2 6.4 6.6 6.8 2 1.2.7 J 1 : 0 i. i.3 0 1 2 3 i: 3 _3 _2 _1 0 1 2 3 1 >: 1. <: 1 %. matrix divide = 12
+/ +/ 1 2 3 6 13
1.3 y t k y = C 00 + C 01 t + + c 0k t k + ɛ S 0 S 1 S k S 1 S 2 S k+1 S 2 S 3 S k+2...... S K S K+1 S 2k c 0 c 1 c 2 c K T 1 T 2 T 3 T k * 10 Example 2 y = 4.89554 + 1.13288x 0.141456x 2 + 0.00438882x 3 3 poly0 DN12 4.89554 1.13288 _0.141456 0.00438882 12 10 8 6 4 2 0-2 0 5 10 15 20 3 ( ) ( )4 4 linefit_poly0 DN12 *10 14
10 3 poly0 10?. 20 _11.2 17.6678 _3.75583 0.227855 f x = 11.2 + 17.6678x 3.75583x 2 + 0.227855x 3 X,.a = 10?. 20 X 0 1 2 3 20 10 X 1 (>:i. # a),.a 1 6 2 3 3 19 4 15 5 10 6 14 7 0 8 7 9 12 10 17 1 3 0 1 2 3 0 (1 (>:i. # a)ˆ/ i. >: 3 1 1 1 1 1 2 4 8 1 3 9 27 1 4 16 64 1 5 25 125 1 6 36 216 1 7 49 343 1 8 64 512 1 9 81 729 1 10 100 1000 15
AIC AIC k AIC 3 1.3.1 J ˆ?,?. # n $ 1.4 AR 1.4.1 AR (AR Auto Regressive Model) M AR Autoregressive Model S 11 S 12 S 1k S 12 S 22 S 2k S 13 S 23 S 3k...... S 1K S 2K S kk b 1 b 2 b 3 b K = T 1 T 2 T 3 T k Yull Walker Burk Householder *11 (DN12) 3 *11 3 Burk Yull Walker householder 16
1. x t 1, x t 2, x t 3,, x t n (X) 2. (y) y 3. X X NR, y 3 y, x t 0, x t 1, x t 2, x t 3 1 3.4 2 5.5 3 7.7 4 10.2 5 8.5 6 7.6 7 10.1 8 11.5 9 6.9 10 2... 25 10.7 26 10.8 (3}.DN12),.}:."1 >3<\DN12 Y x t 1 x t 2 x t 3 3.4 0 0 0 cut 5.5 3.4 0 0 cut 7.7 5.5 3.4 0 cut 10.2 7.7 5.5 3.4 8.5 10.2 7.7 5.5 7.6 8.5 10.2 7.7 10.1 7.6 8.5 10.2 11.5 10.1 7.6 8.5 6.9 11.5 10.1 7.6 2 6.9 11.5 10.1 10.7 10.2 8.8 7.2 10.8 10.7 10.2 8.8 cut ini f ix(\) X Box( ) in f ix(\) Open( ) x t 1, x t 2, x t 3 Rotate(. 1) Y 3 1. "1>3<\ >:i.5 3 2 1 4 3 2 5 4 3 (%.) 17
1.4.2 AIC Mll = n k n (logq k/(n k)) AIC = 2 MLL + 2 k = (n k) (logq k /(n k)) + 2 k AR AIC k 4 exam_ar0 DN12 +----------+-------------+-----------+ mean=6.904 corr=0.546834 AIC=41.4362 +----------+-------------+-----------+ 3 exam_ar0 DN12 +----------+-------------+-----------+ mean=6.904 corr=0.536786 AIC=41.5126 +----------+-------------+-----------+ 5 exam_ar0 DN12 +----------+------------+-----------+ mean=6.904 corr=0.56474 AIC=41.7886 +----------+------------+-----------+ 4 ar0 DN12 0.973496 _0.548562 0.232687 0.0641936 4 0.973496t 1 0.548562t 2 + 0.232687t 3 + 0.0641936t 4 4 linefit_ar DN12 pd eps \temp\ar_0.eps 18
12 native estim 10 8 6 4 2 0 0 2 4 6 8 10 12 1.5 reg0 reg exam ad reg0 reg exam ad n estim reg0 line f it reg0 estim reg0 n linefit reg0 n poly0 poly exam linefit poly0 ar0 exam ar0 4 poly exam n 4 linefit poly0 n m ar0 n m exam ar0 n 19
2 2.1 2.1.1 ]a=: >: i.10 1 2 3 4 5 6 7 8 9 10 Script Example 1 n n X i am=: +/ % # i=1 am2=: # % +/ 5.5 am a n Π n i=1 X i gm=: # %: */ n ni=1 1 X i hm=: am &.(%"_) gm a 4.52873 0 */ 0 hm a 3.41417 20
common mean m x m 1 + + x m n n cm=:[: {.(am,gm)ˆ:_ cm a 5.00257 ( ) 1 f 1 ni=1 f x i n 2.1.2 - km km km km km km km hm 30 40 60 40 3 % +/ 1r30 1r40 1r60 40 3 1 4+3+2 120 = 360 9 2.2 2.2.1 (( y t y t 1 ) 4 1) 100 2.2.2 qtr_grow 1 1.02 NB. 1 1.02 8.24322 1.0 1.02 8.2% 2.2.3 Script qtr_grow=: 3 : 100 * <: ˆ&4 %/. y NB. e.g. u 1 1.02 NB. exchange rate of growth from quarter to year 21
(rotate.) %/. 1 1.02 1.02 4 ˆ4 1 (<: 1 ) <: ˆ&4 %/. 1 1.02 0.0824322 100 ( times) 100 * <: ˆ&4 %/. 1 1.02 8.24322 2.2.4 % n xt+n x t 1 100 Script grow_ave=: 4 : 100 * <: (x %: %/. y) %/. 504827 539160 NB. GDP 1995/2000 1.06801 5 5 %: %/. 504827 539160 NB. GDP 1995/2000 1.01325 22
100 100* 5 %: %/. 504827 539160 NB. GDP 1995/2000 101.325 Working Example 5 grow_ave 504827 539160 NB. GDP 1995/2000 1.32463 1.32% 10 grow_ave 469567 539160 NB. GDP 1990/2000 1.39161 (1.39%) 2.2.5 - DN20 1.25 1.4 1.07143 1.06667 1.125 gm DN20 1.17608 Working Example DN21 1.042 1.125 1.063 1.073 1.134 year 01 02 03 04 05 % 4.2 12.5 6.3 7.3 13.4 gm DN21 1.08681 NB. 8.67% per year J 23
i. 0 >: 1 * % *: square 2 %: square root n {. From _ in f inity child "1 Rank 1 ˆ: power &. under y. Rotate / Insert ( ) <: Decriment 1 & ( ) r 1r3 1 3 24
2.3 2000 100 = y t y 00 3 Laspiress pt q 0 p0 q 0 Parshe pt q t p0 q t Fischer pt q 0 p0 q 0 pt q t p0 q 0 DN22 A B +------+-----+ 75 50 60 70 NB. 85 48 55 72 NB. 1 95 46 52 75 NB. 2 105 44 48 80 NB. 3 +------+-----+ DN22=:(75 85 95 105,. 50 48 46 44) ;60 55 52 48,. 70 72 75 80 ( ) lsp_chain0 DN22 1 1.01887 1.05535 1.08302 NB. 1 2 3 par_chain0 DN22 25
1 1.01515 1.04025 1.04444 fis_chain0 DN22 1 1.01701 1.04777 1.06356 Working Example 1 DN23=:13 170 15 155,52 49 55 41,13 89 14 85,51 57 53 53,41 55 35 57,:45 33 48 31 DN23=: ( ;1 0 1 0 )<;.1 DN23 DN23 2000 2005 +------+------+ 13 170 15 155 NB. Cabbages 52 49 55 41 NB. Spinachs 13 89 14 85 NB. Napa 51 57 53 53 NB. Leek 41 55 35 57 NB. Lettuce 45 33 48 31 NB. Broccoli +------+------+ 2 2 14 lsp DN23 Las: Par: Fis: 103.654 103.209 103.431 2005 2.3.1 Script lsp=: 3 : 0 NB. Calc Laspi Parshe Fischer P0 Q0 P1 Q1 =: {;("2),. : L:0 y las=.(+/ P1 * Q0)% +/ P0 * Q0 par=.(+/p1 * Q1)% +/ P0 * Q1 fis=.%:(las * par) Las: Par: Fis:,: 9j3 ":100 *(las, par,fis ) 26
) Laspeyres Étinne Laspeyres(1834-1913) Lass-pey-ress 26 1871 Wikipedia J n { take Las character 12j5 ": f ormat (12 5 ) 2.4 x x 1 (x x) 2 n x cm y (%) DN24 5.6 30 5.8 26 6 33 6.2 31 6.4 33 6.4 35 6.4 37 6.6 36 6.8 33 stand DN24 _1.7781 _0.843274 _1.22628 _2.10819 _0.674453 0.105409 _0.122628 _0.527046 0.429198 0.105409 0.429198 0.737865 0.429198 1.37032 0.981023 1.05409 1.53285 0.105409 27
1 0 key fruits water plot : stand DN24 pd eps \temp\kanaya_02.eps -1-2 0 1 2 3 4 5 6 7 8 4 ( ) 2.4.1 Script stand=: dev % "1 sd dev=: -"1(+/ % #) sd=. %:@var var=: # % ([:+/[: *: dev) J *: square 2 (ˆ2 ) %: square root 2 ( ) -"1(+/%#) x x [: Cap 28
2.4.2 x x -(+/%#) ( ) 2.5 e log ln J ˆ. ˆ.100 ˆ 4.60517 ˆ 4.60517 100 e y = x logx = y 1x1 2.71828 J 1x1 e ˆ. DN24 1.72277 3.4012 1.75786 3.2581 1.79176 3.49651 1.82455 3.43399 1.8563 3.49651 1.8563 3.55535 1.8563 3.61092 1.88707 3.58352 1.91692 3.49651 29
J ˆ. naturel log 10 ˆ. 100 is 2 ˆ e n ˆ 4.60517 is 100 1x1 e e 2x1 1x2 2.6 x = x x = x1 2 + x2 2 + + x2 n ( ) 2 1 (Euclid) norm=: [: %: [: +/ *: 2.6.1 The length(or norm)of v is the nonnegative scalar v x = v v = x1 2 + x2 2 +... + x2 n v 2 = v v x y = (x 1 y 1 ) 2 + (x 2 y 2 ) 2 (x n y n ) 2 PQ Example euc_norm 1 _2 2 0 0.333333 _0.666667 0.666667 0 v = 1 2 2 0 v = (1 2 2 2 2 2 0) = 9u = 1 v v 2.6.2 Script norm=: [: %: [: +/ *: euc_norm=: 3 : y. % norm y. Coffee Brake John Napier(1550 1617) 1614 e 30
2.7 2.7.1 x = 1 n n i=1 X i (x x) 2 ss=: [: +/ (*:@dev) mean=: +/ % # NB. am dev=: - mean var=:ss%# Variance (x x) 2 n Standard deviation (x x) 2 sd x n sd=: %:&(ss%#) vr=: sd % mean cov=: # % ([: +/ [: */"1 dev) 1 (x x)(y ȳ) N dev=: -"1(+/ % #) cov=:# % ([: +/ [:*/ "1 (-"1(+/ % #))) ( ) 31
2.7.2 V(X 1 + X 2 + + X n ) = V = n V(X i + Cov(X i, X j ) i=1 S xx S xy S xz S yx S yy S yz S zx S zy S zz i j x y x1 y1 x2 y2 x3 y3 x4 C = y4 t XX N x 1 x y 1 ȳ = 1 x 1 x x 2 x... x n x x 2 x y 2 ȳ N y 1 ȳ y 2 ȳ... y n ȳ x 3 x y 3 ȳ x n x y n ȳ ( ) ( 1 (xn x)(x n x) (x n x)(y n ȳ) x = N (y n ȳ)(x n x) (y n ȳ)(y n ȳ) ) 32
2.7.3 Worked Example DN25 NB. (1) 145 30 60 70 145 35 70 75 150 35 65 80 150 40 70 70 155 40 75 75 155 45 70 80 160 40 80 85 160 50 70 90 165 40 65 90 165 45 75 85 dev2 DN25 NB. (2) _10 _10 _10 _10 _10 _5 0 _5 _5 _5 _5 0 _5 0 0 _10 0 0 5 _5 0 5 0 0 5 0 10 5 5 10 0 10 10 0 _5 10 10 5 5 5 33
( : dev2 DN25) +/. * (dev2 DN25) NB. (3) 500 275 175 425 275 300 150 250 175 150 300 100 425 250 100 500 n n = 10 (( : dev2 DN25) +/. * (dev2 DN25)) % # DN25 NB. (4) 50 27.5 17.5 42.5 27.5 30 15 25 17.5 15 30 10 42.5 25 10 50 vartable DN25 50 27.5 17.5 42.5 27.5 30 15 25 17.5 15 30 10 42.5 25 10 50 2.7.4 Script vartable=:# % :@dev2 +/.* dev2 J +/. * : Transpose 2.8 R = 1 r xy r xz r yx 1 r yz r zx r zy 1 X, Y ρ XY = Cov(X, Y) V(X) V(Y) ρ xy = (x, y) x y 2.8.1 Working Example cortable DN25 34
1 0.710047 0.451848 0.85 0.710047 1 0.5 0.645497 0.451848 0.5 1 0.258199 0.85 0.645497 0.258199 1 2.8.2 Script cortable=: 3 : 0 ss=. [: +/ [: *: dev2 sd=. %:&(ss%#) stand=. dev2@] %"1 sd@] cortable=. #@] % ( :@stand@] +/. * stand@]) cortable y ) J @ Atop -"1(+/%#) -"1 is hook 2.9,,,, r i j R = [ r i j ] R R 1 r i j r i j.o = ri j r ii r j j X 1 X 2 R = [ r i j ] R R 1 r i j r i j.o = ri j r ii r j j 35
5 DN26 1 198005 2854 17662 2 254020 5880 24208 3 290237 11995 30538 4 348839 19605 36568 5 452356 23490 48113 cortable DN26 1 0.967571 0.997407 0.967571 1 0.976533 0.997407 0.976533 1 pcor_table cortable DN26 1 _0.414878 0.965817 _0.414878 1 0.631025 0.965817 0.631025 1 0.63 2.10 stand stand n mean dev var sd cov 36
vartable vartable n cortable cortable n pcor table pcor table cortable n am am2 gm hm am i.10 gm i.10 hm i.10 qtr grow grow ave qtr grow 1 1.02 5 grow ave 504827 539160 lsp lsp n lsp chain0 par chain0 fis chain0 norm norm n 37
A ] a=. i. 5 4 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 from ( 0 1 3 _1 { a 0 1 3 _1 {a NB. rank is default 0 1 2 3 4 5 6 7 12 13 14 15 16 17 18 19 from ( 0 2 _1 {"1 a ( 1) 0 2 _1 {"1 a 0 2 3 4 6 7 8 10 11 12 14 15 16 18 19 take ({.) 2{."1 a 0 1 4 5 8 9 12 13 16 17 38
B EXCEL EXCEL Libre CALC biff-8 EXCEL2003 EXCEL2007 biff-12 biff XML ) EXCEL2003 (2007 OK,2010?) tara Libre-Office () B.1 tara Net J Run/Package Manager tables/excel,tables/tara DL Net J CDROM,Net J602 copy B.2 EXCEL ( 0-99999 B.3 tara require files B.3.1 tara.ijs tara.ijs j602/addons/tables/tara/tara.ijs tutorial (tara.ijt) addons/tables/tara/tara.ijt. dir =. /data/sna/esri/principal/2010/ a=.readexcel dir, shouhi_test.xls tara Open BOX ;("1) 9}. 2 4 {"1 a Sheet. 39
Sheet1 readexcel dir, test_calc.xls +--+--+--+-+ 1 2 3 +--+--+--+-+ 2 3 4 +--+--+--+-+ 4 2 6 +--+--+--+-+ 45 65 34 +--+--+--+-+ B.3.2 *12 EXCEL bi=. conew biffbook writenumber bi 0 0 ;i. 10 10 writenumber bi 0 0 ;a1 NB. (example) a1=.? 10 10 $ 100 save bi /temp/testtara.xls underbar 2 ) a1=. i.4 5 1 a1 writexlsheets /temp/tararest.xls *12 tara jmacros.xls J602 csv 40
B.4 CSV B.4.1 CSV save CSV Comma Separated Values EXCEL ( 0-99999 copy EXCEL csv Example CSV 1994 http://www5.cao.go.jp/keizai3/getsurei.html index shouhi test.csv save B.4.2 CSV J require files csv dir=: c:/data/sna/esri/principal/2010/ ] a=. readcsv dir, shouhi_test.csv ] a=. ".@> readcsv dir, shouhi_test.csv ] a=. ;("1) ".(L:0) a 41
References 2000 J 1996 [] 1985 Miscellance J602 is download available (No charge) http://www.jsoftware.com Scripts are accessible http://japla.sakura.ne.jp 42