# NONPAR v2.1

#

# Nonparametric estimation of the density and empirical distribution for a time series

#

# 2 data series: (1) simulated data and (2) the quarterly growth rate of the US Gross National Product

#

# For the second series we specifically have the Real Gross National Product, USA, quarterly 1947-Q1 to 2006-Q2, FRED database series GNPC96

# Federal Reserve Bank of Saint Louis online database http://research.stlouisfed.org/fred2/, in billions of chained 2000 dollars

#

# Initialize

clear

ginit

clc

ranseed

# Simulate a nonlinear autoregressive model, pp. 16-17 of reference book

N = 1400.

zt = zeros(N,1)

et = unirnd(-1.,1.,N,1)

for i = 2 to N begin

zt(i) = 2*zt(i-1)/(1+0.8*zt(i-1)^2) + et(i)

end

# Estimate and plot the density using the second reference bandwidth

[grid,f,h] = density(zt,1,-3.5,+3.5,201,"hopt2",0,1)

# Now plot the simulated series and the density in a single window

gaddwin()

up = gaddview(0.,0.5,1.,0.5)

gxaxis("linear",0,N,10)

gxgrid("major")

gygrid("major")

gtitle("Simulated Nonlinear Autoregressive Time Series")

gplot(zt)

down = gaddview(0.,0.,1.,0.5)

gxgrid("major")

gygrid("major")

gcolor("red")

gtitle(["Kernel Density estimated with bandwidth h = ",ntoa(h)])

gplot(grid,f)

# Now estimate and plot the empirical distribution function and compare it with the standard normal

#

# Note that the distribution of zt is not normal but a Kolmogorov-Smirnov test cannot reject the null

# hypothesis of normality (the p-value is, however, obtained under independence). Try a larger sample value,

# for example N = 1500 to see how the result changes

#

std_zt = (zt-colmean(zt))/colstd(zt)

[grid,F] = distribution(std_zt,-3.5,+3.5,101,0)

gaddwin()

gtitle("EDF for Simulated Nonlinear Autoregressive Time Series (black) and the N(0,1) distribution (red)")

gxgrid("major")

gygrid("major")

gyaxis("linear",0,1,10)

gcolor({"black","red"})

gplot(grid,[F,cnormal(grid)])

print "P-value of Kolmogorov-Smirnov test (simulated) = ",kolsmi(std_zt,function cnormal)

print

# Read US data and divide by 1e+4

N = 238.

usgnp = read("data\usgnp4706.prn","double",N,1)/1e+4

# Convert to quarterly growth rate

gusgnp = diff(log(usgnp))*100

# Estimate and plot the density using the second reference bandwidth

[grid,f,h] = density(gusgnp,1,-4.5,+4.5,101,"hopt2",0,1)

# Now plot the series and the density in a single window

gaddwin()

up = gaddview(0.,0.5,1.,0.5)

gxaxis("linear",0,N+2,10)

gxgrid("major")

gygrid("major")

gtitle("Quarterly Growth Rate (in %) of the Real US Gross National Product, 1947-2006")

gplot(gusgnp)

down = gaddview(0.,0.,1.,0.5)

gxgrid("major")

gygrid("major")

gcolor("red")

gtitle(["Kernel Density estimated with bandwidth h = ",ntoa(h)])

gplot(grid,f)

# Now estimate and plot the empirical distribution function and compare it with the standard normal

std_gusgnp = (gusgnp-colmean(gusgnp))/colstd(gusgnp)

[grid,F] = distribution(std_gusgnp,-4.5,+4.5,101,0)

gaddwin()

gtitle("EDF of the Quarterly Growth of the US GNP (black) and the N(0,1) distribution (red)")

gxgrid("major")

gygrid("major")

gyaxis("linear",0,1,10)

gcolor({"black","red"})

gplot(grid,[F,cnormal(grid)])

print "P-value of Kolmogorov-Smirnov test (US GNP)= ",kolsmi(std_gusgnp,function cnormal)

print