-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathSimulationQuantileReg.R
143 lines (118 loc) · 3.43 KB
/
SimulationQuantileReg.R
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
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
Mu <- 5.00
N <- 100
epsilon <- rnorm(N,0,1)
Y_t <- Mu + epsilon
hist(Y_t)
plot(Y_t,type='p')
postBeta <- function(Beta,Y,N,p){
post <- N*log(p) + N*log(1-p) - sum((abs(Y-Beta)+((2*p-1)*(Y-Beta)))/2)
return(post)
}
updateBeta <- function(Beta,Y,N,p,Sm,Sigma,t,BetaMean){
accept <- NULL
count <- 0
c <- .8
ggamma <- 1/(t^c)
proposed <- Beta + sqrt(Sm*Sigma)*rnorm(1,0,1)
probab <- min(1,exp(postBeta(proposed,Y,N,p)-postBeta(Beta,Y,N,p)))
if(runif(1)<probab){
accept <- proposed
count <- 1
} else {
accept <- Beta
}
lsm <- log(Sm)+ggamma*(probab - 0.234)
Sm <- exp(lsm)
Sigma <- Sigma + ggamma*(((Beta - BetaMean)^2)-Sigma)
BetaMean <- BetaMean + ggamma*(Beta-BetaMean)
return(list(accept,count,Sm,Sigma,BetaMean))
}
##########################################################################
Niter <- 20000
N <- 100
Beta.out1 <- array(NA, dim = Niter)
Count1 <- array(0, dim = Niter)
Sm1 <- array(NA, dim = Niter)
Sigma1 <- array(NA, dim = Niter)
BetaMean1 <- array(NA, dim = Niter)
Beta.out1[1] <- .5
Sm1[1] <- 2.4^2
Sigma1[1] <- 2
BetaMean1[1] <- 3
Beta.out2 <- array(NA, dim = Niter)
Count2 <- array(0, dim = Niter)
Sm2 <- array(NA, dim = Niter)
Sigma2 <- array(NA, dim = Niter)
BetaMean2 <- array(NA, dim = Niter)
Beta.out2[1] <- .5
Sm2[1] <- 2.4^2
Sigma2[1] <- 2
BetaMean2[1] <- 3
Beta.out3 <- array(NA, dim = Niter)
Count3 <- array(0, dim = Niter)
Sm3 <- array(NA, dim = Niter)
Sigma3 <- array(NA, dim = Niter)
BetaMean3 <- array(NA, dim = Niter)
Beta.out3[1] <- .5
Sm3[1] <- 2.4^2
Sigma3[1] <- 2
BetaMean3[1] <- 3
Beta.out4 <- array(NA, dim = Niter)
Count4 <- array(0, dim = Niter)
Sm4 <- array(NA, dim = Niter)
Sigma4 <- array(NA, dim = Niter)
BetaMean4 <- array(NA, dim = Niter)
Beta.out4[1] <- .5
Sm4[1] <- 2.4^2
Sigma4[1] <- 2
BetaMean4[1] <- 3
t1 <- Sys.time()
for(i in 2:Niter){
Beta1 <- updateBeta(Beta.out1[i-1],Y_t,N,.05,Sm1[i-1],Sigma1[i-1],i,BetaMean1[i-1])
Beta.out1[i] <- Beta1[[1]]
Count1[i] <- Beta1[[2]]
Sm1[i] <- Beta1[[3]]
Sigma1[i] <- Beta1[[4]]
BetaMean1[i] <- Beta1[[5]]
Beta2 <- updateBeta(Beta.out2[i-1],Y_t,N,.25,Sm2[i-1],Sigma2[i-1],i,BetaMean2[i-1])
Beta.out2[i] <- Beta2[[1]]
Count2[i] <- Beta2[[2]]
Sm2[i] <- Beta2[[3]]
Sigma2[i] <- Beta2[[4]]
BetaMean2[i] <- Beta2[[5]]
Beta3 <- updateBeta(Beta.out3[i-1],Y_t,N,.75,Sm3[i-1],Sigma3[i-1],i,BetaMean3[i-1])
Beta.out3[i] <- Beta3[[1]]
Count3[i] <- Beta3[[2]]
Sm3[i] <- Beta3[[3]]
Sigma3[i] <- Beta3[[4]]
BetaMean3[i] <- Beta3[[5]]
Beta4 <- updateBeta(Beta.out4[i-1],Y_t,N,.95,Sm4[i-1],Sigma4[i-1],i,BetaMean4[i-1])
Beta.out4[i] <- Beta4[[1]]
Count4[i] <- Beta4[[2]]
Sm4[i] <- Beta4[[3]]
Sigma4[i] <- Beta4[[4]]
BetaMean4[i] <- Beta4[[5]]
print(i)
}
t2 <- Sys.time()
(t2-t1)
plot(Beta.out1[1001:20000],type='l')
plot(Beta.out2[1001:20000],type='l')
plot(Beta.out3[1001:20000],type='l')
plot(Beta.out4[1001:20000],type='l')
hist(Beta.out1[1001:20000])
hist(Beta.out2[1001:20000])
hist(Beta.out3[1001:20000])
hist(Beta.out4[1001:20000])
Beta1 <- mean(Beta.out1[1001:20000])
Beta2 <- mean(Beta.out2[1001:20000])
Beta3 <- mean(Beta.out3[1001:20000])
Beta4 <- mean(Beta.out4[1001:20000])
plot(Y_t,type='p')
abline(h=Beta1,col='red')
abline(h=Beta2,col='red')
abline(h=Beta3,col='red')
abline(h=Beta4,col='red')
meansVec <- c(Beta1,Beta2,Beta3,Beta4)
quants <- quantile(Y_t,probs=c(.05,.25,.75,.95))
cbind(meansVec,quants)