# Graphs

- Shade U.S. recessions as identified by NBER (ShadeUSRecessions.prg)
- Plotting joint confidence ellipse (joint.prg)
- Plotting empirical distribution functions (edfplot.prg)

## Shade U.S. recessions as identified by NBER

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

## Plotting joint confidence ellipse

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

## Plotting empirical distribution functions

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)