Using BasisMatrices.jl
to reproduce the figures in
Miranda and Fackler, Applied computational economics and finance, Chapter 6.
In [1]:
using BasisMatrices
using Plots
using DataFrames
In [2]:
plotlyjs()
Out[2]:
In [3]:
a, b = 0, 1
x = collect(linspace(a, b, 1001));
In [4]:
n = 10
basis = Basis(ChebParams(n, a, b))
Out[4]:
In [5]:
ys = []
m = 9
for i in 1:m
c = zeros(n)
c[i] = 1
y = funeval(c, basis, x)
push!(ys, y)
end
ylims = (-1.1, 1.1)
plot(x, ys, layout=m, ylims=ylims, leg=false)
Out[5]:
In [6]:
k = 1
m = 9
breaks = m - (k-1)
basis = Basis(SplineParams(breaks, a, b, k))
Out[6]:
In [7]:
ys = []
for i in 1:m
c = zeros(m)
c[i] = 1
y = funeval(c, basis, x)
push!(ys, y)
end
ylims = (-0.1, 1.1)
plot(x, ys, layout=m, ylims=ylims, leg=false)
Out[7]:
In [8]:
k = 3
m = 9
breaks = m - (k-1)
basis = Basis(SplineParams(breaks, a, b, k))
Out[8]:
In [9]:
ys = []
for i in 1:m
c = zeros(m)
c[i] = 1
y = funeval(c, basis, x)
push!(ys, y)
end
ylims = (-0.1, 1.1)
plot(x, ys, layout=m, ylims=ylims, leg=false)
Out[9]:
Consider the univariate function $f(x) = \exp(-\alpha x)$ on $[-1, 1]$ as in Section 6.7.
In [10]:
f(x, alpha) = exp(-alpha * x)
a, b = -1, 1;
Chebychev approximation:
In [11]:
alpha = 1
n = 10
basis = Basis(ChebParams(n, a, b))
S, (xgrid,) = nodes(basis)
Phi = BasisMatrix(basis, Expanded(), S, 0)
y = f.(xgrid, alpha)
c = Phi.vals[1] \ y
x = collect(linspace(a, b, 1001))
interp = funeval(c, basis, x);
Approximation error:
In [12]:
ylims = (-1e-9, 1e-9)
plot(x, f.(x, alpha)-interp, ylims=ylims, yformatter = :scientific, leg=false)
Out[12]:
Consider $f_1(x) = \exp(-x)$, $f_2(x) = (1 + 25x^2)^{-1}$, and $f_3(x) = \lvert x\rvert^{0.5}$ as in Section 6.6.
In [13]:
f_1(x) = exp(-x)
f_2(x) = 1./(1 + 25x.^2)
f_3(x) = sqrt(abs(x));
In [14]:
a, b = -5, 5
x = collect(linspace(a, b, 1001))
degrees = [10, 20, 30]
errors = Array{Float64}(3, 3, 3);
In [15]:
methods = [1, 2] # Linear Spline, Cubic Spline
ks = [1, 3]
for (method, k) in zip(methods, ks)
for (j, degree) in enumerate(degrees)
breaks = degree - (k-1)
basis = Basis(SplineParams(breaks, a, b, k))
S, (xgrid,) = nodes(basis)
Phi = BasisMatrix(basis, Expanded(), S, 0)
for (i, func) in enumerate([f_1, f_2, f_3])
y = func.(xgrid)
c = Phi.vals[1] \ y
interp = funeval(c, basis, x)
errors[i, j, method] = maximum(abs, interp - func.(x))
end
end
end
In [16]:
method = 3 # Chebychev Polynomial
for (j, degree) in enumerate(degrees)
basis = Basis(ChebParams(degree, a, b))
S, (xgrid,) = nodes(basis)
Phi = BasisMatrix(basis, Expanded(), S, 0)
for (i, func) in enumerate([f_1, f_2, f_3])
y = func.(xgrid)
c = Phi.vals[1] \ y
interp = funeval(c, basis, x)
errors[i, j, method] = maximum(abs, interp - func.(x))
end
end
In [17]:
col_names = map(Symbol, ["Degree", "Linear Spline", "Cubic Spline", "Chebychev Polynomial"]);
$\exp(-x)$:
In [18]:
df = DataFrame(vcat(Any[degrees], Any[errors[1, :, method] for method in 1:3]), col_names)
Out[18]:
$(1 + 25x^2)^{-1}$:
In [19]:
df = DataFrame(vcat(Any[degrees], Any[errors[2, :, method] for method in 1:3]), col_names)
Out[19]:
$\lvert x\rvert^{0.5}$:
In [20]:
df = DataFrame(vcat(Any[degrees], Any[errors[3, :, method] for method in 1:3]), col_names)
Out[20]:
In [21]:
degree = 7
a, b = -1, 1
x = collect(linspace(a, b, 1001));
Chebychev polynomial:
In [22]:
basis_cheb = Basis(ChebParams(degree, a, b))
S, (xgrid_cheb,) = nodes(basis_cheb)
Phi_cheb = BasisMatrix(basis_cheb, Expanded(), S, 0);
Cubic spline:
In [23]:
k = 3
breaks = degree - (k-1)
basis_cub = Basis(SplineParams(breaks, a, b, k))
S, (xgrid_cub,) = nodes(basis_cub)
Phi_cub = BasisMatrix(basis_cub, Expanded(), S, 0);
Linear spline:
In [24]:
k = 1
breaks = degree - (k-1)
basis_lin = Basis(SplineParams(breaks, a, b, k))
S, (xgrid_lin,) = nodes(basis_lin)
Phi_lin = BasisMatrix(basis_lin, Expanded(), S, 0);
In [25]:
bases = [basis_cheb, basis_cub, basis_lin]
xgrids = [xgrid_cheb, xgrid_cub, xgrid_lin]
Phis = [Phi_cheb, Phi_cub, Phi_lin];
Approximation of $\exp(-x)$:
In [26]:
func = f_1
ys = [func.(x)]
for (basis, xgrid, Phi) in zip(bases, xgrids, Phis)
y = func.(xgrid)
c = Phi.vals[1] \ y
interp = funeval(c, basis, x)
push!(ys, interp)
end
ylims = (0, 3)
plot(x, ys, layout=length(ys), xlims=(a, b), ylims=ylims, leg=false)
Out[26]:
Approximation of $(1 + 25x^2)^{-1}$:
In [27]:
func = f_2
ys = [func.(x)]
for (basis, xgrid, Phi) in zip(bases, xgrids, Phis)
y = func.(xgrid)
c = Phi.vals[1] \ y
interp = funeval(c, basis, x)
push!(ys, interp)
end
ylims = (-0.2, 1.2)
plot(x, ys, layout=length(ys), xlims=(a, b), ylims=ylims, leg=false)
Out[27]:
Approximation of $\lvert x\rvert^{0.5}$:
In [28]:
func = f_3
ys = [func.(x)]
for (basis, xgrid, Phi) in zip(bases, xgrids, Phis)
y = func.(xgrid)
c = Phi.vals[1] \ y
interp = funeval(c, basis, x)
push!(ys, interp)
end
ylims = (0, 1)
plot(x, ys, layout=length(ys), xlims=(a, b), ylims=ylims, leg=false)
Out[28]:
In [ ]: