*! version 0.1.0, Brent McSharry, 12july2013 * Risk adjusted exponentially weighted moving average chart program define raewma, sortpreserve version 10.1 syntax varlist(min=2 max=2 numeric) [if] [in] , Predicted(varname numeric) STartest(real) /* */ [SMooth(real 0.01) Alpha1(real 0.05) Alpha2(real 0.01) XLABEL(passthru) YTITLE(passthru) YLABEL(passthru) LEGEND(passthru) RESolution(real 0.0000125)] tempvar ewma_o ewma_p ewma_se cl1_ub cl1_lb cl2_ub cl2_lb cl touse mark `touse' `if' `in' markout `touse' `varlist' `predicted' local observed:word 1 of `varlist' local timevar:word 2 of `varlist' if (`smooth'<=0 | `smooth'>1) { di as error "smooth option must be (0,1]" return 125 } if (`resolution'<0 | `resolution'>0.01) { di as error "resolution option must be [0,0.01]" return 125 } if (`startest'<0 | `startest'>1) { di as error "smooth option must be [0,1]" return 125 } qui count if `touse' local nobs `r(N)' if (`nobs'==0) { di as error "nobs observed" return 2000 } gsort -`touse' `timevar' mata: getRaewma("`predicted'","`observed'",`nobs',`smooth',`startest',"`ewma_o'","`ewma_p'","`ewma_se,'", `resolution') qui { local z=invnormal(1-(`alpha2'/2)) gen `cl' = `z'*`ewma_se' gen `cl2_ub' = `ewma_p' + `cl' gen `cl2_lb' = `ewma_p' - `cl' local z=invnormal(1-(`alpha1'/2)) replace `cl' = `z'*`ewma_se' gen `cl1_ub' = `ewma_p' + `cl' gen `cl1_lb' = `ewma_p' - `cl' } list `cl2_ub' if `touse' & _n>5380 twoway (rarea `cl2_lb' `cl2_ub' `timevar' if `touse', color(gs11)) /* */ (rarea `cl1_lb' `cl1_ub' `timevar' if `touse', color(gs13)) /* */ (line `ewma_o' `timevar' if `touse', pstyle(p1)) /* */ , `xlabel' `ytitle' `ylabel' legend(rows(1) order(3 2 1) label(1 "`=100-`alpha2'*100'% CI") label(2 "`=100-`alpha1'*100'% CI") label(3 "Observed") `legend') end mata: void getRaewma(string scalar expectVarname, string scalar obsVarname, real scalar usedNobs, real scalar smooth, real scalar startVal, string scalar emwaObsName, string scalar emwaExpectName, string scalar emwaSeName, real scalar resolution) { real colvector ewmaExp, ewmaObs, ewmaSe real rowvector smoothPow,usedPow real vector expectVar, colIndices real scalar oneMinusSm, j, k, len, currExp, currObs, expectJ, cumSeCalc, oIndx, eIndx, endPow, startExpect oneMinusSm = 1 - smooth oIndx = st_varindex(obsVarname) eIndx = st_varindex(expectVarname) smoothPow = getSmoothPow(smooth, usedNobs, resolution) endPow = cols(smoothPow) startExpect = 1 len = st_nobs() ewmaExp = J(len,1,.) ewmaObs = J(len,1,.) ewmaSe = J(len,1,.) expectVar = J(usedNobs,1,.) currExp = currObs = startVal for (j=1;j<=usedNobs;j++) { currObs = ewmaObs[j,1] = _st_data(j,oIndx)*smooth + currObs *oneMinusSm expectJ = _st_data(j,eIndx) currExp = ewmaExp[j,1] = expectJ*smooth + currExp *oneMinusSm cumSeCalc = expectVar[j] = expectJ*(1-expectJ) if (j<=endPow) { usedPow = smoothPow[|1,(endPow+1-j) \ 1,endPow|] } else { ++startExpect } ewmaSe[j,1] = smooth*sqrt( usedPow * expectVar[|startExpect,1 \ j,1|]) } (void) setbreakintr(0) colIndices = st_addvar("double", (emwaExpectName, emwaObsName, emwaSeName), 1) st_store(.,colIndices[1],ewmaExp) st_store(.,colIndices[2],ewmaObs) st_store(.,colIndices[3],ewmaSe) (void) setbreakintr(1) } //build lookup for power calculations (which tend to be computationally expensive) real rowvector function getSmoothPow(real scalar smooth, real scalar length, real scalar resolution) { real scalar i, oneMinusSm, oneMinusSmSq, currPow, nocols real rowvector returnVar oneMinusSm = (1-smooth) oneMinusSmSq = oneMinusSm^2 nocols = ceil(ln(1-(1-resolution)^2*(1-oneMinusSmSq^(length+1)))/ln(oneMinusSmSq))-1 if (nocols > length) { nocols = length } returnVar = J(1,nocols,.) currPow = returnVar[1,nocols] = 1 for (i=nocols-1;i>0;i--) { returnVar[1,i] = currPow = oneMinusSmSq * currPow } return(returnVar) } end