-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathgsreg_example-1.jl
88 lines (63 loc) · 2.32 KB
/
gsreg_example-1.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
# gsres function example
# cd("C:\\Users\\millsjf\\OneDrive - University of Cincinnati\\Class\\9011 2022")
## data here:
using Distributions, Plots, StatsPlots
n = 1000
x = randn(n)
y = 1.0 .+ 0.5.*x .+ 0.4.*randn(n)
plot(x,y,st=:scatter)
## AR(1) DGP
phi = 1.001
z = zeros(n+20)
for t = 2:(n+20)
z[t] = 0.0 + phi*z[t-1] + randn(1)[1]
end
y = z[21:(n+20)]
plot(y)
yt = y[2:end]
yt1 = y[1:end-1]
###############################################
include("gsreg.jl")
## gsreg(y::Array{Float64},X::Array{Float64}; M=10000::Int64, burnin = Int(floor(M/10.0))::Int64,tau=[1.0]::Array{Float64},iB0=[0.0001]::Array{Float64},b0=[0.0]::Array{Float64},d0=0.0001::Float64, a0=3.0::Float64)
# uninformative prior
b = [0.0; 0.0] # prior coeff. means
iB = inv([0.0001 0.0; 0.0 0.0001])
X = [ones(n-1) yt1]
bdraws,s2draws = gsreg(yt,X)
plot(bdraws[:,2], st=:density, fill=true, alpha=0.5, title = "phi posterior")
mean(bdraws[:,2])
std(bdraws[:,2])
quantile(bdraws[:,2],[0.025,0.975])
plot(s2draws[:,1], st=:density, fill=true, alpha=0.5, title = "sigma^2 posterior" )
mean(s2draws[:,1])
std(s2draws[:,1])
quantile(s2draws[:,1],[0.025,0.975])
# Informative prior (still uninformative for variance parameter)
b = [0.0; 1.0] # prior coeff. means
iB = inv([1000.0 0.0; 0.0 0.001]) # prior cov matrix
bdraws,s2draws = gsreg(y,X, b0=b, iB0=iB)
plot!(bdraws[:,2],st=:density,linecolor=:red,label="informative")
mean(bdraws[:,2])
std(bdraws[:,2])
quantile(bdraws[:,2],[0.025,0.975])
plot(s2draws, st=:density, fill=true, alpha=0.5, title = "σ squared posterior", label= false )
z = randn(1000,2).*0.3 .+ 0.5
plot(z[:,1], z[:,2], st = :scatter) #, xlims = (-0.2,1.2), ylims = (-0.2,1.2), label=false)
#vline!([0.0 1.0], label=false)
#hline!([0.0 1.0], label=false)
# correlated variables
z3 = 0.5*z[:,1] .+ randn(1000).*0.2
cor(z[:,1], z3)
plot(z[:,1], z3, st = :scatter, xlabel = "z1", ylabel="z3")
# AR(1) model estimation
# uninformative prior
b = [0.0; 0.0] # prior coeff. means
iB = inv([0.0001 0.0; 0.0 0.0001])
yt = y[2:n]
yt1 = y[1:(n-1)]
X = [ones(n-1) yt1]
bdraws,s2draws = gsreg(yt,X)
plot(bdraws[:,2], st=:density, fill=true, alpha=0.5, title = "phi posterior", label= "uninformative" )
mean(bdraws[:,2])
std(bdraws[:,2])
quantile(bdraws[:,2],[0.025,0.975])