Graphs
The following program shows some of the command line options to
customize the look of your graph. It calls a subroutine to add shades
for the periods identified as U.S. recessions by the NBER. The source
can be found at
US Business Cycle Expansions and Contractions.
' shade U.S. recessions identified by NBER
' for version 4
' last revised 10/27/2000 h
'include subroutine file
include sub_ShadeUSRecessions.prg
'change path to program path
%path = @runpath
cd "{%path}"
' create workfile
workfile tmp1 q 1948:1 1998:4
' fetch data from database
fetch(d=data_svar) rgnp unemp
' plot data
group g1 dlog(rgnp) unemp
freeze(graph1) g1.line
' set legend text
graph1.setelem(1) legend(Real GNP growth)
graph1.setelem(2) legend(Unemployment rate)
' dual scaling
graph1.setelem(2) axis(right)
' set frame size
graph1.options size(8,2)
' draw zero line
graph1.draw(dashline, left, rgb(172,172,172)) 0
' shade recessions
call shadeUSRecessions(graph1, 4)
show graph1
The subroutine file sub_shadeUSRecessions.prg looks as follows:
' shade recessions identified by NBER
' 9/19/2000 h
'
' g1: name of graph object to add shades
' f : scalar for frequency (f=4 for quarterly and f=12 for monthly)
subroutine shadeUSRecessions(graph g1, scalar f)
if (f = 12) then
g1.draw(shade,bottom) 1857:06 1858:12
g1.draw(shade,bottom) 1860:10 1861:06
g1.draw(shade,bottom) 1865:04 1867:12
g1.draw(shade,bottom) 1869:06 1870:12
g1.draw(shade,bottom) 1873:10 1879:03
g1.draw(shade,bottom) 1882:03 1885:05
g1.draw(shade,bottom) 1887:03 1888:04
g1.draw(shade,bottom) 1890:07 1891:05
g1.draw(shade,bottom) 1893:01 1894:06
g1.draw(shade,bottom) 1895:12 1897:06
g1.draw(shade,bottom) 1899:06 1900:12
g1.draw(shade,bottom) 1902:09 1904:08
g1.draw(shade,bottom) 1907:05 1908:06
g1.draw(shade,bottom) 1910:01 1912:01
g1.draw(shade,bottom) 1913:01 1914:12
g1.draw(shade,bottom) 1918:08 1919:03
g1.draw(shade,bottom) 1920:01 1921:07
g1.draw(shade,bottom) 1923:05 1924:07
g1.draw(shade,bottom) 1926:10 1927:11
g1.draw(shade,bottom) 1929:08 1933:03
g1.draw(shade,bottom) 1937:05 1938:06
g1.draw(shade,bottom) 1945:02 1945:10
g1.draw(shade,bottom) 1948:11 1949:10
g1.draw(shade,bottom) 1953:07 1954:05
g1.draw(shade,bottom) 1957:08 1958:04
g1.draw(shade,bottom) 1960:04 1961:02
g1.draw(shade,bottom) 1969:12 1970:11
g1.draw(shade,bottom) 1973:11 1975:03
g1.draw(shade,bottom) 1980:01 1980:07
g1.draw(shade,bottom) 1981:07 1982:11
g1.draw(shade,bottom) 1990:07 1991:03
else
if (f = 4) then
g1.draw(shade,bottom) 1857:2 1858:4
g1.draw(shade,bottom) 1860:4 1861:2
g1.draw(shade,bottom) 1865:2 1867:4
g1.draw(shade,bottom) 1869:2 1870:4
g1.draw(shade,bottom) 1873:4 1879:1
g1.draw(shade,bottom) 1882:1 1885:2
g1.draw(shade,bottom) 1887:1 1888:2
g1.draw(shade,bottom) 1890:3 1891:2
g1.draw(shade,bottom) 1893:1 1894:2
g1.draw(shade,bottom) 1895:4 1897:2
g1.draw(shade,bottom) 1899:2 1900:4
g1.draw(shade,bottom) 1902:3 1904:3
g1.draw(shade,bottom) 1907:2 1908:2
g1.draw(shade,bottom) 1910:1 1912:1
g1.draw(shade,bottom) 1913:1 1914:4
g1.draw(shade,bottom) 1918:3 1919:1
g1.draw(shade,bottom) 1920:1 1921:3
g1.draw(shade,bottom) 1923:2 1924:3
g1.draw(shade,bottom) 1926:4 1927:4
g1.draw(shade,bottom) 1929:3 1933:1
g1.draw(shade,bottom) 1937:2 1938:2
g1.draw(shade,bottom) 1945:1 1945:4
g1.draw(shade,bottom) 1948:4 1949:4
g1.draw(shade,bottom) 1953:3 1954:2
g1.draw(shade,bottom) 1957:3 1958:2
g1.draw(shade,bottom) 1960:2 1961:1
g1.draw(shade,bottom) 1969:4 1970:4
g1.draw(shade,bottom) 1973:4 1975:1
g1.draw(shade,bottom) 1980:1 1980:3
g1.draw(shade,bottom) 1981:3 1982:4
g1.draw(shade,bottom) 1990:3 1991:1
else
statusline currently only monthly or quarterly frequency supported!
endif
endif
endsub
The following program plots the joint confidence region (an ellipse) of two
parameters estimated by OLS. The program estimates a multiple linear regression,
extracts two parameters into a vector and its covariance matrix into a sym
matrix and calls a subroutine to plot an ellipse.
' plot joint confidence region of two OLS parameters
' replicates Greene (1993,p.191)
' for version 4
' 12/20/2000 h
' include subroutine
include drawEllipse.prg
' load workfile
%data = @runpath + "greene6_2"
load "{%data}"
' estimate equation by OLS
equation eq1.ls y x1 x2 x3 x4 x5
' get c(2) and c(3)
vector(2) beta
beta(1) = eq1.c(2)
beta(2) = eq1.c(3)
' get covariance martrix of c(2) & c(3)
sym(2,2) binv
binv(1,1) = eq1.@cov(2,2)
binv(1,2) = eq1.@cov(2,3)
binv(2,2) = eq1.@cov(3,3)
' and invert
binv = @inverse(binv)
' call the subroutine to plot the joint confidence ellipse
!fcv95 = @qfdist(0.95, 2, eq1.@regobs-eq1.@ncoef)
call drawEllipse(beta, binv, 2*!fcv95, "gra1") ' 3rd arg is critical value
' edit legend and title to graph
gra1.setelem(1) legend(b_2)
gra1.setelem(2) legend(b_3)
gra1.addtext(t) Joint vs Individual 0.95 Confidence Region
' add individual 0.95 region as shades
!b_low = eq1.c(2)-2*eq1.@stderrs(2)
!b_upp = eq1.c(2)+2*eq1.@stderrs(2)
gra1.draw(shade,bottom) !b_low !b_upp
gra1.draw(dashline,bottom) !b_low
gra1.draw(dashline,bottom) !b_upp
!b_low = eq1.c(3)-2*eq1.@stderrs(3)
!b_upp = eq1.c(3)+2*eq1.@stderrs(3)
gra1.draw(shade,left) !b_low !b_upp
gra1.draw(dashline,left) !b_low
gra1.draw(dashline,left) !b_upp
show gra1
The subroutine to draw the ellipse looks as follows:
' subroutine to plot an ellipse represented as a quadratic form
' (x-x0)'B(x-x0) = c
'
' where the inputs are
'
' x0: 2x1 vector
' B: 2x2 symmetric matrix
' c: scalar
' %gname: name for graph object
'
' 1/22/99 h.
' revised 12/19/2000 for version 4 (h)
subroutine drawEllipse(vector x0, sym B, scalar c, string %gname)
' input error check
!rows = @rows(B)
!cols = @columns(B)
if (@rows(x0)<>2 or !rows<>2 or !cols<>2) then
statusline input arguments must have dimension 2
return
endif
' get eigenvalues/eigenvectors
matrix ve = @eigenvectors(B)
vector va = @eigenvalues(B)
' construct the rotation matrix
!ang = @atan(ve(2,1)/ve(1,1)) ' angle
matrix(2,2) rot
rot.fill(b=r) @cos(!ang),-@sin(!ang),@sin(!ang),@cos(!ang)
' draw the unrotated ellipse
vector(101) theta
!pi = @acos(-1)
for !i = 1 to 101
theta(!i) = (!i-1)*2*!pi/100
next
!a = @sqrt(c/va(1)) ' width of ellipse
!b = @sqrt(c/va(2))
matrix(2,101) xy
for !i=1 to 101
xy(1,!i) = !a*@cos(theta(!i))
xy(2,!i) = !b*@sin(theta(!i))
next
' rotate
matrix rx = rot*xy
' offset
matrix ones = @filledmatrix(1,101,1)
rx = rx + x0*ones
rx = @transpose(rx)
' plot
freeze({%gname}) rx.scat
' connect scatter without symbols
{%gname}.setelem(1) lpat(solid) symbol(none)
endsub
Version 4 includes built-in tests for goodness-of-fit based on empirical
distribution functions. The program edfplot.prg illustrates how to
obtain graphical output to visually assess the goodness-of-fit. In particular,
it illustrates how to overlay the kernel density plot with the theoretical
density under test.
'plots for goodness-of-fit tests
'for version 4
'12/18/2000 h.
'--------------------------------------------------------------------------------------
'subroutines
'--------------------------------------------------------------------------------------
'evaluate weibull density at xt and store in fx
subroutine eval_weibull(series xt, series fx, scalar m, scalar s, scalar a)
series zt = (xt - m)/s
fx = (a/s) * zt^(a-1) * exp( -(zt^a) )
endsub
'evaluate Gaussian density at xt and store in fx
subroutine eval_norm(series xt, series fx, scalar m, scalar s)
series zt = (xt - m)/s
!pi = @acos(-1)
fx = (1/@sqrt(2*!pi)/s) * exp( -(0.5*zt^2) )
endsub
'--------------------------------------------------------------------------------------
'main program
'--------------------------------------------------------------------------------------
'overlay theoretical and empirical distribution
workfile edfplot u 1 500
'set random number generator
rndseed(type=mt) 1234567
'simulate data
series x = @rchisq(5)
'kernel density estimate
!ngrid = 150 'number of points to evaluate
do x.kdensity(!ngrid, b=2.3, o=kdxy)
'make sure workfile is large enough
if (@obsrange < !ngrid) then
expand 1 !ngrid
endif
'convert stored estimates into series
smpl 1 !ngrid
series xt
series ft_kdensity
group group1 xt ft_kdensity
mtos(kdxy, group1)
'evaluate theoretical distribution at xt
series ft_dist
call eval_norm(xt, ft_dist, @mean(x), @stdev(x))
graph g1.xyline xt ft_dist ft_kdensity
g1.option linepat ' need to set linepat
g1.setelem(1) legend()
g1.setelem(2) legend(normal)
g1.setelem(3) legend(kernel density)
g1.setelem(1) lcolor(blue) lpat(dash1)
g1.setelem(2) lcolor(red) lpat(solid)
g1.legend columns(1) -inbox position(2.3,0)
'qq-plot
freeze(g2) x.qqplot(n)
g2.setelem(1) legend()
g2.addtext(0.1,0.1) QQ-plot
'combine two graphs
graph graph1.merge g2 g1
graph1.align(2,0,0)
show graph1
'display goodness-of-fit tests
show x.edftest(type=normal,showopts)