ECON546 lab8 Bayesian Estimation

1. Data

import data from web source

# import data from web source
cons.df<- read.table("http://web.uvic.ca/~dgiles/blog/consump.dat",header=T)

data, demeaned, scaled

consump<- (cons.df$CONS-mean(cons.df$CONS))/1000
income<- (cons.df$Y-mean(cons.df$Y)) /1000
year<- cons.df$YEAR

2. MLE estimator

ols and mle, non constant

# ols and mle, non constant
mle<- lm(consump~income -1)
summary(mle)

Call:
lm(formula = consump ~ income - 1)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.043304 -0.002994  0.001686  0.008586  0.047164 

Coefficients:
       Estimate Std. Error t value Pr(>|t|)    
income  0.89848    0.00581   154.7   <2e-16 ***
---
Signif. codes:  0 ?**?0.001 ?*?0.01 ??0.05 ??0.1 ??1

Residual standard error: 0.01869 on 35 degrees of freedom
Multiple R-squared:  0.9985,    Adjusted R-squared:  0.9985 
F-statistic: 2.392e+04 on 1 and 35 DF,  p-value: < 2.2e-16

prepare data for later computation

c2<-consump^2
y2<-income^2
cy<-consump*income

3. Prior

Prior for \(\pi\) marginal propencity of consumption: mean and variance

\[\mu_{\pi} = 0.8\]

\[ \sigma^2_{\pi} = 0.001\]

pm<- 0.8
pv<- 0.001
n<- length(consump)

deduce prior parameters for \(\beta\)

Joint prior distribution

\[ p(\pi) \sim \mathrm{Beta}(\pi| \alpha, \beta) \]

\[ p(\sigma) \propto 1/\sigma \]

We know that

\[\bar{\pi}=\alpha/(\alpha + \beta)=0.8 \]

and

\[ \sigma^2_{\pi} = \frac{\alpha\beta}{(\alpha + \beta)^{2}(\alpha + \beta + 1)} = 0.0001 \]

See page 2

 
alpha<- (1/pv)*pm*(pm*(1-pm)-pv)
gamma<- (1/pv)*(1-pm)*(pm*(1-pm)-pv)

Prior pdf for Marginal propencity of consumption

curve(dbeta( x, shape1 =  alpha, shape2 = gamma), 
      from = 0.64, to = 0.92 ,col="red",xlab= "Pi",
     ylab="Density", main="Prior pdf for MPC" ,lwd=3)

      
#       
#  curve(dgamma(x,shape=s+a,rate=n+b),col=”red”,xlab=”LAMBDA”,
# +   ylab=”DENSITY”,lwd=3,from=3,to=25)
# > curve(dgamma(x,shape=a,rate=b),col=”blue”,lwd=3,add=TRUE)
# > curve(dgamma(x,shape=s+1,rate=n),col=”green”,lwd=3,add=TRUE)
# > legend(“topright”,c(“PRIOR”,”LIKELIHOOD”,”POSTERIOR”),
# +  col=c(“blue”,”green”,”red”),lty=1,lwd=3)      
    

4. Joint Posterior

On page 4 we get joint posterior, then what we can do is to get marginal Posteriors for \(b\) and \(\sigma\)

\[p(\pi, \sigma | c)\]

And then we can get the Bayes’ point estimators by using some loss functions

\(c, c^2, y^2 \ \text and \ cy\) are data

4.1 First integrate through sigma to get a marginal posterior for \(\pi\)

\(\sigma\) can take value from 0 to \(\infty\)

Then we get

\[p(\pi |c)\]

make a function which can be fed in R using native integrate function from page 5/6

integrand<- function(b){(b^(alpha-1))*((1-b)^(gamma-1))/(sum(c2)+b^2*sum(y2)-
                                                       2*b*sum(cy))^(n/2)}

to get normalizing constant \(k\) page 6

(norm<- 1/integrate(integrand, lower = 0, upper = 1)$value)

[1] 5089.043

4.2 Posterior mean for MPC

under quadratic loss, function for mean of the posterior

posterior mean, integrate through beta from 0 to 1, page 7

page 7

\[\int_{0}^{1} \pi p(\pi |c) d \pi \]

See there is one more b in the front of integrand

integrand<- function(b) {b*(b^(alpha-1))*((1-b)^(gamma-1))/(sum(c2)+b^2*sum(y2)-
                                                              2*b*sum(cy))^(n/2)}
(post_mean<- norm*integrate(integrand, lower = 0, upper = 1)$value)

[1] 0.8931085

The “posterior uncertainty” about MPC \(\pi\) can be measured, for instance, by computing

Notice (b-post_mean)^2 in the integrand, same pdf and same normalized constant.

#  posterior variance,  page 7
integrand<- function(b) {(b-post_mean)^2*(b^(alpha-1))*
    ((1-b)^(gamma-1))/(sum(c2)+b^2*sum(y2)-2*b*sum(cy))^(n/2)}
(post_var<- norm*integrate(integrand, lower = 0, upper = 1)$value)

[1] 1.115461e-05

4.4 Prior

# generate prior densicy
b<- seq(0.68 , 0.93, 0.00001)
prior<- dbeta(b,alpha, gamma)

4.5 Likelihood

\[ p(c | \pi, \sigma) \propto \left((n-1)\sigma^2 + (\pi - \pi_{mle})^2 \sum y^2 \right)^{\frac{n-1}{2}} \]

#  mle estimators for slope/mpc: beta and standard diviation: sigma 
sig<- summary(mle)$sigma
bmle<- summary(mle)$coef[1]
# function for integration of likelihood
integrand<- function(b) {((n-1)*sig^2+(b-bmle)^2*sum(y2))^(-(n-1)/2)} 
(norml<-1/integrate(integrand, lower = 0, upper = 1)$value)

[1] 2.272411e-32

# likelihood/ density function 
mlf<- norml*((n-1)*sig^2+(b-bmle)^2*sum(y2))^(-(n-1)/2)

4.6 marginalize the Joint Posterior for \(\sigma\) by integrating out \(\pi\).

to get

\[ p(\sigma | c)\]

See page 8 and 9, but we do not do it in this lab

4.7 Plot Prior, Likelihood, and Posterior

posterior of MPC , see page 6

# posterior of MPC , see page 6
posterior<- norm*((b)^(alpha-1))*((1-b)^(gamma-1))/(sum(c2)+b^2*sum(y2)
-2*b*sum(cy))^(n/2)

Mode of posterior of MPC

which( posterior ==   max(posterior)) 

[1] 21354

posterior[which( posterior ==   max(posterior)) ]

[1] 69.435

# calculate posterior mode
o<- order(posterior)
temp<- b[o]
(post_mode<- temp[25001])

[1] 0.89353

Comparision

plot(b, posterior, type="l", lwd=2, col="blue",xlab="Beta", ylab="Densities")
lines(b, prior, lwd=2,col="red")

lines(b, mlf, lwd=2, col="black")
#abline(v = post_mean, col = "blue", lwd =2 )
title(main="Consumption Function Example")

mtext("(Prior Mean = 0.8 ; Prior Variance = 0.001)")
legend(0.7,50,lty=c(1,1,1), lwd=c(2,2,2), col=c("red","black", "blue"), c("Prior","Likelihood","Posterior"))

Different point estimators

bmle # MLE of Beta

[1] 0.8984838

post_mean # Posterior Mean for Beta

[1] 0.8931085

post_var # Posterior Variance for Beta

[1] 1.115461e-05

post_mode # Posterior Mode for Beta

[1] 0.89353