// Nov 13 2011 11:38:45 // Based on bstrpr.do, attempt to generalise as an ado // Represent the average time-dependent transition pattern underlying a // set of sequences capture program drop trprgr // VARLIST is the set of consecutive variables describing the sequences (must be reshape-able) // Options // - required // NSTATES is number of states (numbered consecutively from 1) // - optional // Probabilities on the diagonal (p(i,i)) are plotted with range FLOOR to 1 // Probabilities off the diagonal (p(i,j), i!=j) are plotted with range 0 to CEILING // GMAX is the highest count in any state in any period // MOVINGAVERAGE is the half-window width for the smoothing in the transition matrix caclulation // TEXTSIZE is a textsizestyle option that controlls the textsize program define trprgr syntax varlist(min=2) [if] [in], /* NStates(int) */ [FLoor(real 0.9) CEIling(real 0.1) GMax(integer 400) MOVingaverage(int 3) TEXtsize(string) *] marksample touse local gmaxxs `gmax' local floorxs `floor' local ceilingxs `ceiling' // allow textsizestyle options if ("`textsize'"!="") { local textsize ",size(`textsize')" } else { local textsize ",size(`textsize')" } preserve keep if `touse' tempname id gen `id' = _n keep `id' `varlist' local state : word 1 of `varlist' local statelab: value label `state' local state = regexr("`state'","[0-9]+","") local seql : word count `varlist' tempname timevar tempname trpr // Store the smoothed transition matrix di "Creating transition probability matrix" qui maketrpr `varlist', mat(`trpr') ma(`movingaverage') // Reshape long to calculate the // chronograms di "Creating chronogram data" qui reshape long `state', i(`id') j(`timevar') qui tab `state', gen(m) qui sum `state' local nstates `r(max)' collapse (sum) m1-m`nstates' , by(`timevar') label variable `timevar' "Time" di "Drawing chronograms" // Get intervals for count scale; principle copied from `gmax'units in hmap.ado local step 1 while `gmax' == round(`gmax',`step') { local step = `step'*10 } /* i is such that `gmax'==round(`gmax',i) is not true but `gmax'==round(`gmax',i/10) is */ /* e.g. `gmax' is 1234500, i is 1000. */ local step = `step'/10 if (`gmax'!=round(`gmax',1)) { local step 1 while `gmax' == round(`gmax',`step') { local step = `step'/10 } local step = `step'*10 } forvalues i = 1/`nstates' { if ("`statelab'"!="") { local sl : label `statelab' `i' label variable m`i' "`sl'" } else { label variable m`i' "State `i'" } qui su m`i' if (r(max)>`gmaxxs') { local gmaxxs `r(max)' noisily display in blue "Warning: state `i' exceeds GMAX: `r(max)' > `gmax'" } local cl : word `i' of navy maroon forest_green dkorange teal cranberry lavender khaki sienna emidblue emerald brown erose gold bluishgray tempname tmpgraph`i' twoway area m`i' `timevar', nodraw pstyle(p`i') yscale(range(0 `gmax')) xtitle(`textsize') ytitle(`textsize') ylabel(#5) /* ylabel(0(`step')`gmax') */ `options' name(`tmpgraph`i'', replace) /* graph save `tmpgraph`i'', replace */ } local gcom "" //clear qui drop _all qui svmat `trpr' qui gen period = 1 + int((_n-1)/`nstates') qui gen state = 1 + mod(_n-1,`nstates') di "Drawing transition time-series" forvalues origin = 1/`nstates' { local gcom "`gcom' `tmpgraph`origin''" if ("`statelab'"!="") { local yt : label `statelab' `origin' } else { local yt "Origin `origin'" } forvalues dest = 1/`nstates' { if ("`statelab'"!="") { local xt: label `statelab' `dest' } else { local xt "Destination `dest'" } tempname graphof_`origin'_`dest' local gcom "`gcom' `graphof_`origin'_`dest''" if (`origin'==`dest') { qui su `trpr'`dest' if state==`origin' if (r(min)<`floorxs') { local floorxs `r(min)' noisily display in blue "Warning: Min t-rate(`origin',`dest') is less than floor: `r(min)' < `floor'" } line `trpr'`dest' period if state==`origin', lwidth(*2) nodraw yscale(range(`floor' 1)) ylabel(`floor'(0.025)1.0) ytitle("`yt'"`textsize') xtitle("`xt'"`textsize') `options' name(`graphof_`origin'_`dest'', replace) } if (`origin'!=`dest') { qui su `trpr'`dest' if state==`origin' if (r(max)>`ceilingxs') { local ceilingxs `r(max)' noisily display in blue "Warning: Max t-rate(`origin',`dest') is greater than ceiling: `r(max)' > `ceiling'" } line `trpr'`dest' period if state==`origin', nodraw yscale(range(0 `ceiling')) ylabel(0.0(0.025)`ceiling') ytitle("`yt'"`textsize') xtitle("`xt'"`textsize') `options' name(`graphof_`origin'_`dest'', replace) } /* graph save `graphof_`origin'_`dest'', replace */ } } if (`gmaxxs' > `gmax') { di "GMAX exceeded: `gmaxxs'" } if (`floorxs' < `floor') { di "FLOOR exceeded: `floorxs'" } if (`ceilingxs' != `ceiling') { di "CEILING exceeded: `ceilingxs'" } di "Combining graphs" if ("`options'" != "") { local options ", `options'" } graph combine `gcom' `options' restore end