NB. script for numeric recipe for econometrics NB. chapter 2 3 4 5 NB. written by Masato Shimura 14/Mar/2006 NB. jcd02773@nifty.ne.jp NB. Ver 1.0 NB. ----------------UTIL sel=: 2 : 'm. {"1 n.' NB. over and by is famous script in J over=:({.;}.)@":@, by=:' '&;@,.@[,.] byover=: 3 : ' ( i. # y ) by (i. {: $ y ) over y ' byover2=: 3 : ' y by y over y ' NB. arrange yourself NB. =============================== NB. chapter 3 NB. ================================================================ NB. D. Durbin-Watson Test NB. ================================================================ dw=: 4 : 0 NB. usage: x dw y ( not append 1: to y ) a=. x %. 1 ,. y znew=. (1 ,. y ) +/ . * a b=. +/ (tmp=:x - znew)^2 NB. u^2 u1=. 2 <\ tmp u2=. (}: "1 >u1) - {: "1 >u1 NB. ut - ut-1 u3=. +/ u2 ^2 ans=. u3 % b ) NB. =============================== NB. chapter 4 NB. ============================ NB. calculate mav mav12 NB. ============================ mean=: +/%# dev=: -(+/%#) dev2=: -"1 (+/%#) mav=: +/\ % [ NB. Season adjustment using 12 month mav12=: [: 2&mav 12&mav mavm=: +/\("1) % ("1) [ NB. mav for matrix : usagge: 12 mavm data(matrix yoko) mav_c=:2&mav@mav ` mav @. (2&|@[) NB. ================================= NB. apply for vertical matrix vmav=: +/\ % "1 [ vmav12=: [: 2&vmav 12&vmav NB. =================================== mav_w3=: 3 : 0 wt=: 1 2 3 2 1 % 9 M0=:> +/ L:0(5<\ y ) * L:0 wt ) NB. complement lack of w_mav3 NB. y is yel mav data cpl_mavw3=: 3 : 0 Y=: mav_w3 y wt_0=: 3 7 10 7 % 27 wt_1=: 2 3 4 % 9 HEAD0=: +/ (|. wt_1)* 3 {. y HEAD1=: +/(|. wt_0) *4 {. y TAIL99=:+/ wt_1 * _3{. y TAIL98=:+/ wt_0 * _4{. y TMP=: (HEAD0,HEAD1),Y,(TAIL98,TAIL99) ) NB. ==================================== NB. plot combination data and mav12 NB. Usage: plotcomb y NB. ==================================== plot_mav12=: 3 : 0 NB. plot mav 12 and actual data NB. Ussage: plotcomb_12 y tmp0=:(7#6),7 + i. # tmp1=.mav12 y tmp0=.tmp0 ,7 #({: tmp0) NB. x axis tmp2=:( 7# {. tmp1), tmp1, 7 # {: tmp1 NB. y axis tmp3=.tmp0;tmp2 pd 'reset' pd y pd 'color red' pd tmp3 pd 'show' ) mav12cpl=: 3 : 0 NB. mav12 and hokan by centered simple mav NB. amari sinnrai sinai houga yoi NB. Usage: mav12cpl y BAI=: 0.5,(11#1),0.5 NB. 0.5 1 1 ..1 0.5 BAI2=: 6}. }: |. L:0 <\ BAI NB. /. .1 1 0.5 TMP0=:12 %~ L:0 +/ L:0 BAI * (L:0) 13<\y FST=: 12{. y NB. early times FST1=:(0.5 -~ L:0 # L:0 FST0)%~ L:0 +/ L:0 BAI2 * L:0 FST0=: 6}. 12{. <\FST NEAR=: |. _12{.y NB. recent times FST3=:(0.5 -~ L:0 # L:0 FST2)%~ L:0 +/ L:0 (|. L:0 BAI2) * L:0 FST2=: 6}. 12{. <\ NEAR (>FST1),(>TMP0),|. > FST3 ) mav_w13=: 3 : 0 NB. WEG=: _11 0 9 16 21 24 25 24 21 16 9 0 -11 WEG=: _11 0 9 16 21 24 25 24 21 16 9 0 _11 TMP0=: 143 %~ L:0 +/ L:0 WEG * (L:0) 13<\y WEGN=: |. (L:0)_6{. }: <\ WEG WEGF=: _6{.}: <\ WEG NB. 5}. |. L:0 L:0 <\ WEG FST=: 12{. y NB. early times FST1=:(+/ L:0 WEGF)%~ L:0 +/ L:0 WEGF * L:0 FST0=: 6{. 6}. <\FST NEAR=: |. _12{.y NB. recent times FST3=:(+/ L:0 WEGN)%~ L:0 +/ L:0 WEGN * L:0 FST2=: 6{. 6}. <\ NEAR (>FST1),(>TMP0),|. > FST3 ) plot_w13=: 3 : 0 NB. compare mav12cpl /fine!! pd 'reset' pd y pd 'color green ' pd mav12cpl y pd 'color red' pd mav_w13 y pd 'show' ) plot_mav12cpl=: 3 : 0 pd 'reset' pd y pd 'color red' pd mav12cpl y pd 'show' ) NB. ================================= NB. mav any dimension NB. ussage: x plotcomb y (x dimension y data) NB. ================================== p_sub0=: 4 : 0 NB. make mav anytime odd or even mav=. +/\ % [ mav_even=.[: 2&mav mav select. tmp0 =.2 | x case. 0 = tmp0 do. goto_odd. case. 1 = tmp0 do. goto_even. end. label_odd. tmp10=.x mav y tmp2=.<. -: x return. label_even. tmp10=. x mav_even y tmp2=.>. -: x return. ) p_sub1=: 4 : 0 NB. make plot data by mav x p_sub0 y tmp30=:(tmp2 # tmp2) , tmp2 +i. # tmp10 tmp3=.tmp30,tmp2 #({: tmp30) NB. x axis tmp4=:( tmp2 # {. tmp10), tmp10,tmp20=.tmp2 # {: tmp10 NB. y axis tmp5=.tmp3;tmp4 ) plotcomb=: 4 : 0 NB. usage: x plotcomb y NB. x jisuu NB. y data pd 'reset' pd x p_sub1 y pd y pd 'color red' pd tmp5 pd 'show' ) plot_w3=: 3 : 0 tmp=. mav_w3 y NB. tmp1=:(2# (2{tmp1))(0 1)} tmp1=: i. # y NB. tmp1=: tmp1 ,2# {: tmp1=: _2}. tmp1 NB. compaund X NB. tmp2=: ((2# {. tmp),tmp,2#{:tmp) NB. expand Y pd 'reset' pd 'color red' pd y pd 'color blue' pd cpl_mavw3 y NB. tmp1;tmp2 pd 'color green' NB. pd spencer y NB. may display spencer pd 'show' ) plot_ws=: 3 : 0 NB. plot w3 and spencer tmp=. mav_w3 y NB. tmp1=:(2# (2{tmp1))(0 1)} tmp1=: i. # y NB. tmp1=: tmp1 ,2# {: tmp1=: _2}. tmp1 NB. compaund X NB. tmp2=: ((2# {. tmp),tmp,2#{:tmp) NB. expand Y pd 'reset' pd 'color red' pd y pd 'color blue' pd cpl_mavw3 y NB. tmp1;tmp2 pd 'color green' pd spencer y NB. may display spencer pd 'show' ) NB. ---spencer---------------------------------- spencer2=: 3 : 0 NB. Spencer 15 weighted mav WEIGHT=. _3 _6 _5 3 21 46 67 74 67 46 21 3 _5 _6 _3 % 320 y,. ; WEIGHT +/ . * (L:0) 15<\ ( (7 # {.),] ,7 #{:) y ) NB. C.Reiter Fractals Visuarization and J (2nd Edition 2000) wts=: (|.,}.) 74 67 46 21 3 _5 _6 _3%320 locspen=: (+/ . *)&wts spencer=: 15"_ locspen\ (7: # {.),] ,7:#{: plot_spencer=: 3 : 0 NB. Usage: plot_spencer y NB. plot y & spencer y pd 'reset' pd 'type line' pd 'itemcolor 255 0 0' NB. Blue pd y NB. walk pd 'type line' pd 'itemcolor 0 0 255' NB. Akane pd spencer y NB. spencer walk pd 'show' ) NB. ========================== NB. henderson moving average NB. this Script is written by J.Takeuchi NB.=========================== hnwgt=: 4 : 0 m=: -: x +3 NB. (n+3)/2 nm=: 315*(((m-1)^2)-*:y )*((m^2)-*:y )*(((m+1)^2)-*:y )*(((3**:m)-16)-11**:y ) hi=:nm%8*m*((m^2)-1)*((4**:m)-1)*((4**:m)-9)*((4**:m)-25) ) henderson=: 4 : 0 NB. use x under kisuu/not guusuu NB. Usage: i.e. 15&henderson y //or 7 henderson y wt=: x hnwgt i:-:<:x x (+/ . *)&wt\((<.@-:x )#{.y ),y ,(<.@-:x )#{:y ) NB. =================== NB. written by M. Shimura 2003/Dec-2004/Mar plot_henderson=: 3 : 0 pd 'reset' pd 'color green' pd y pd 'color blue' NB. pd 'key red=henderson 31' NB. pd 'key blue =Spencer 15' pd spencer y pd 'color red' pd 15&henderson y pd 'show' ) plot_henderson2=:4 : 0 pd 'reset' NB. pd 'pensize 3' pd 'title Henderson 31 and Spencer 15' pd 'color green' pd y pd 'color blue' pd spencer y pd 'color red' pd x henderson y NB. pd 'xlabel 4 5 6 7 8 9 0 1 2 3 ' pd 'keystyle 0' pd 'key y spencer henderson' pd 'show' ) plot_comb_all=: 4 : 0 NB. henderson spencer mav12cpl NB. x is jisuu(henderson)/must kisuu NB. like 31 33 /well as mav12 pd 'reset' pd 'color green' pd y pd 'color blue' pd spencer y pd 'color red' pd x henderson y pd 'color purple' pd mav12cpl y pd 'show' ) NB. make dummy mk_dummy4=: 4 : 0 NB. make dummy for 4 seasons NB. x tag for dummy // select from 0 1 2 3 NB. 0 Spring 1 Summer 2 Autumn 3 Winter // NB. y is data NB. Usage: x mk_dummy4 y Y=. y if. (# $ y )=1 do. Y=. ( i. # y ),. y end. D1=. <|: x =/ 0 1 2 3 TIMES=.>: <.(# Y) % 4 D2=.|. "1 ((# Y),# x ) $ ,> TIMES # D1 ( i. # y ),.D2,.{:"1 Y ) mk_dummy12=: 4 : 0 NB. x is month /start is april NB. x is month to put tag /e.g. 8 1 2 (x mk_dummy12_sub y ),.y ) mk_dummy12_sub=: 4 : 0 TMP0=. ( 3|. >: i. 12) e. x TMP1=. +/ * TMP0 NB. number of raw TMP2=.TMP0,. ( 12, <:TMP1)$ 0 TMP3=.( TMP0 expand =/~ i. TMP1) (i. # y ),. (# y ){. ;(>.12%~ # y ) # < TMP3 ) NB. ===================== NB. exponential smoothing NB. written by Giichiro Suzuki NB.====================== smexp=:4 : 0 s=. {.y=. y while. 1<#y do. s=.s,(x *{:s)+(1-x )*{.y=.}.y end. ) ave_year=: 3 : 0 NB.find month put tag/ rate compare with mean TMP0=._12<\ y NB. classify 12 month TMP1=. (+/ > TMP0) % +/ * > TMP0 NB. +/*> is nr each month TMP2=. ". 6.1 ": 100 * TMP1 % (+/%#) TMP1 NB. (i.12),.(3|. >:i.12),.TMP2 NB. nr month data(%) ) NB. ========================================== NB. complete modified 15/mar/2007 adj_4season_sub=: 3 : 0 NB. find season_index TMP1=.(2}._2}. y) ,. 4&mav_c Y0=: (-4|#y)}. y SEASON0=. 2|. "1 >_4<\ %/"1 TMP1 SEASON1=. +/ SEASON0 % (# TMP1) % 4 SEASON1 * 4 % +/ SEASON1 ) adj_4season=: 3 : 0 NB. using mav_c and reg0 Y0=. (-4|#y)}. y SEASON=. adj_4season_sub Y0 TMP0=. Y0,.;(_4<\Y0) % L:0 SEASON NB. const data f=. reg0 {:"1 TMP0 NB. reg /find ternd TREND=. (1,.>: i. # TMP0) +/ . * f NB. reg_t use >: i. # y TMP0,.TREND,. ;SEASON * (L:0) _4<\ TREND ) plot_4season=: 4 : 0 NB. Usage: e.g. (1994 1) plot_4season a XAXIS=: x mk_4season_sub y pd 'reset' pd 'keypos bo' pd 'keystyle oh' pd 'key data deseason trend fitted' pd XAXIS;|: adj_4season y pd 'show' ) mk_4season_sub=: 4 : 0 NB. mk xaxis 'YY MM'=: x NUM=: # TMP=: _4<\ y YM0=:; 0 0.25 0.5 0.75 + L:0 YY + (L:0) _4<\ 4# i. NUM DROP_TOP=:((>:i.12) e. MM)# 0 0 0 1 1 1 2 2 2 3 3 3 (# y) {.DROP_TOP}. YM0 ) NB. ========12 month========================================= NB. ---------primitive----------------- NB. using whole months to calc seasonal rate NB. update 130 June 2005 /tmp2 add <. NB. slightly modified everywhere 15/Mar/2007 saj=:3 : 0 NB. using whole manthes to calc NB. Season Adjusument using Moving Average NB. Usage: saj y // y is data T0=. (TMP1=:(|(# TMP0)%12 ) , 12) $ TMP0=. mav12 y T1=. TMP1 $ 6}. _6}. y 6 |. T3* 12 % +/ T3=.(+/%#) T1%T0 ) adj_season=: 3 : 0 NB. 12 season data //using centered moving average Y0=. (-(12| # y )) }. y NB. drop out of 4 times from last //if 0 drop nothing S40=.( 6}. _6}. Y0),. mav12 Y0 S41=. 6|."1 >_12<\ %/ "1 S40 NB. _4 non overlap,2|. is rotate S42=.+/ S41 % ({. $ S40) % 12 SEASON=:S42 * 12 %(+/ S42) NB. Sasonal S44=:Y0,. ,> (_12 <\ Y0) % L:0 SEASON NB. Deseasolized: TCe f=. (;{:"1 S44) %. X0=.1 ,. >: i. # S44 NB. f with OLS S44 ,. tmp,.;(_12<\tmp=. X0 +/ . * f)* L:0 SEASON NB. trend ) NB. --------------------------------------- saj2=:4 : 0 NB. last update 15/mar/2006 NB. Season Adjusument using Moving Average NB. x is last month/e.g. 12 NR=: # y if. NR >: 66 do. ANS=. saj_sub0 y NB. whole 5 years elseif. NR >: 54 do. ANS=. saj_sub1 y NB. median elseif. NR < 54 do. ANS=. 12 # 1.0 NB. all 1.0 end. IND=. rotate_month_sub x IND |. ANS ) rotate_month_sub=: 3 : 0 NB. y is last month of data M=. 4 5 6 7 8 9 10 11 12 1 2 3 (11 - M i. y ) - 7 NB. modify 6->7 ) NB. ----------------------------------------------------- saj_check=: 3 : 0 NB. case NR >:60 /drop max min S1=. (_60{. _6}. y ) % S0=. _60{. mav12 y /:~ _12 >\S1 ) saj_sub0=: 3 : 0 NB. case NR >:60 /drop max min S1=. (_60{. _6}. y ) % S0=. _60{. mav12 y S2=.(+/%#) }: }. _12 >\S1 S2 * 12 % +/ S2 ) saj_sub1=: 3 : 0 NB. case. NR>:48 /median S1=. (_36{. _6}. y ) % S0=. _36{. mav12 y S2=. ; 1{. _12 >\S1 NB. median S2 * 12 % +/ S2 ) adj_season2_sub=: 4 : 0 SIND=: (# y ) {. (;(>. (# y ) % 12) # < x saj2 y ) y % SIND NB. SIND is season index ) adj_season2=: 4 : 0 NB. x is last month TMP=. x adj_season2_sub y y ,.TMP,. (1,.i.# y ) +/ . * reg0 TMP ) adj_season_poly=:4 : 0 NB. classical deseason/ 12 month /full length NB. x is (1;12) poly 1 or 2/ last month NB. S1=: (# y ) {. ;(>.(# y )%12) # < (>{: x ) saj2 y NB. full length ( (>{. x ) poly_estim y ) ,. y % S1 ) plot_adj_poly=: 4 : 0 NB. y is x adj_season2 y NB. x is jisuu NB. Usage: x plot_adj_poly 12 adj_season2 y DAT=. (2{."1 y ) ,.{:"1 x poly_estim 1{"1 y pd 'reset' pd |: DAT pd 'show' ) adj_season_reg=:4 : 0 NB. classical deseason/ 12 month /full length NB. x last month e.g. 12 S1=: (# y ) {. ;(>.(# y )%12) # < x saj2 y NB. full length y ,.((reg0 y % S1)&p. i. # y ),. y % S1 ) plot_adj_season=:4 : 0 NB. Usage: e.g. (1994 1) plot_4season a XAXIS=: x mk_season_sub y pd 'reset' pd 'keypos bo' pd 'keystyle oh' pd 'key data deseason trend fitted' pd 'type line' pd XAXIS;}:|: adj_season y NB. pd 'y2axis' NB. pd 'type dot' NB. pd XAXIS;{:"1 adj_season y pd 'show' ) mk_season_sub=:4 : 0 NB. mk xaxis 'YY MM'=. x NUM=. # _12<\ y YM0=.; ((i.12)%12) + L:0 YY + (L:0) _12<\12 # i.>: NUM DROP_TOP=.((>:i.12) e. MM)# >:i.12 (# y) {.DROP_TOP}. YM0 ) NB. plot_adj_season=:4 : 0 NB. DAT=: |: x adj_season2 y NB. pd 'reset' NB. pd 'keypos bottom outside left' NB. pd 'keystyle horizontal' NB. pd 'key data deseason trend' NB. pd 'color green' NB. pd 0{ DAT NB. Y origin NB. pd 'color red' NB. pd 1{ DAT NB. NB. pd 'color blue' NB. pd 2{ DAT NB. trend NB. pd 'show' NB. ) find_rate=: 3 : '(12}. y ) % _12}.y ' NB. very simple! plot_rate=: 3 : 0 RATE=: find_rate y pd 'reset' pd RATE,:(# RATE) # 1 pd 'show' ) NB. ================================ wmav=: 3 : 0 NB. 3ji * 3 ji Weighted MAV NB. Usage: wmav y LS=: (2 # {. y ), y , (2 # {: y ) NB. add 2 + 2 to head and Tailt WT=: 1 2 3 2 1 % 9 NB. Weights T0=: 5 <\ LS T1=: T0 * L:0 WT T2=: > _2}. 2}. +/ L:0 T1 NB. Drop Head and Tail ) sel_3=: 4 : 0 NB. select DATA for plot NB. x is syasyu//0 to 4 NB. y is data/after order_order DAT=:> x {"1 L:0 y SIZE=:$ DAT DAT2=:(i. {. SIZE);(i. {: SIZE);DAT ) NB. =============holt's 2 parameter Smoothing============= holt_loop=:4 : 0 NB. Holt's Method find minimum resideal error NB. Usage: x hw_main y NB. x is Alpha & Beta //y is data//e.g. 0.2 0.6 hw_main NB. use x ex 0.3 0.7 // Less than 1(both) HX0=.10 %~ >. x * 10 NB. take unit 0.1 drop unit 0.01 HX1=.0.05 * >:i. ({. HX0)% 0.05 NB. Alpha to 0 pitch 0.5 HX2=.0.05 * >:i.({: x )% 0.05 HXMAT=. ((# HX2) # HX1),.tmp=. , >(# HX1)# : CTER end. TMP=: HXMAT,.ANS_RES (\:{:"1 TMP){ TMP ) holt_res_sub=:4 : 0 NB. find Residual Standard Error H0=. x holt y H2=. -/"1 H1=:}. 0 3 {"1 H0 RSE=. %: +/ (*: H2 )% # H2 ) holt_main=: 4 : 0 NB. Holt's Method main program NB. Usage: x hw_main y NB. x is Alpha & Beta //y is data//e.g. 0.2 0.6 hw_main NB. use x ex 0.3 0.7 // Less than 1(both) TMP=: x holt_loop y NB. sort NB. HXX=.TMP /: 2{"1 TMP (}: {: TMP) holt y NB. do hw with optimized parameter ) holt=:4 : 0 NB. //calc part NB. Holt's 2 patrameter trend NB. x alpha;beta NB. y data 'AL BE'=. x NB. alpha beta ANS_S=.({. y ) 0} (# y ) # 0 COUNTER=. 1 Y=. 4{. y B=.-: +/ > -/"1 (L:0) |.(L:0) _2<\ Y S=.(AL *{. y )+(1-AL)*( {.ANS_S )+ B ANS_S=. S (1) } ANS_S ANS_B=. B (0 ) }(#y )#0 while. COUNTER < # y do. S=.(AL * COUNTER{y )+(1-AL)*((<:COUNTER) {ANS_S )+(<: COUNTER){ ANS_B ANS_S=. S COUNTER } ANS_S B=. (BE*(S-(<:COUNTER){ANS_S))+(1-BE)* (<: COUNTER){ANS_B ANS_B=. B COUNTER } ANS_B COUNTER=. >: COUNTER end. y ,.ANS_S,. ANS_B,.(}: 0,ANS_S+ANS_B) ) plot_holt=: 4 : 0 NB. same holt plot |: }. 0 3 {"1 x holt y )