### Supporting R code for the linear model elicitation example in Daneshkhah and Oakley (2010)

### Implements the method in
### Kadane, J. B., Dickey, J. M., Winkler, R. L., Smith, W. S. and Peters, S. C. (1980). 
### Interactive elicitation of opinion for a normal linear model. 
### Journal of the American Statistical Association, 75: 845-854.


### Specify design matrix X

x1<-rep(c(20000,70000),4)
x2<-rep(c(250,250,500,500),2)
x3<-c(10,10,10,10,20,20,20,20)
x4<-c(600,200,200,600,600,200,200,600)
x5<-c(4,8,8,4,8,4,4,8)
X<-cbind(rep(1,8),x1,x2,x3,x4,x5)

# Can obtain fractional factorial design using the BHH2 library:
# library(BHH2)
# ffDesMatrix(5,gen=list(c(4,1,2),c(5,3,4)))


### Initial judgements
y5<-matrix(c(385,200,575,260,700,390,765,575),8,1)
y75<-matrix(c(410,235,605,295,735,430,800,620),8,1)
y9<-matrix(c(430,265,630,330,770,470,840,665),8,1)

### Obtain b
b<-solve(t(X)%*%X)%*%t(X)%*%y5

### Compare fitted medians with elicited medians
fitted.medians<-X%*%b
plot(y5,fitted.medians)

### Choose d
a<-(y9-y5)/(y75-y5) #tail ratios. Should inspect for outliers

d.fit<-function(delta,mean.a){
(qt(0.9,delta)/qt(0.75,delta)-mean.a)^2
}

d<-optimise(d.fit,interval=c(0.1,100),mean.a=mean(a))$minimum

### hypothetical data
y0<-matrix(c(417,259,627,306,775,471),6,1)

### Conditional medians
### Element i,j is median for Y_j conditional on y^0_1,...y^0_i
C.C<-matrix(0,6,6)
C.C[1,]<-c(390,205,580,265,705,395)
C.C[2,]<-c(0,210,585,270,710,400)
C.C[3,]<-c(0,0,590,275,715,405)
C.C[4,]<-c(0,0,0,280,720,410)
C.C[5,]<-c(0,0,0,0,730,415)
C.C[6,]<-c(0,0,0,0,0,425)

### Conditional 75th percentiles
### Element i+1 is median for Y_{i+1} conditional on y^0_1,...y^0_i
new.75<-c(0,240,615,310,755,455)

### Estimating R and omega
S.C<-rep(0,6) 
U.i<-matrix(((y75[1]-y5[1])/qt(0.75,d))^2,1,1) # U_1
S.C[1]<-U.i[1,1]

omega.hat<-rep(0,6)

### First estimate of omega
cospread<-(C.C[1,1]-y5[1])*U.i[1,1]/(y0[1]-y5[1])
omega.hat[1]<-U.i[1,1]-cospread


for(i in 1:5){
    
    #S(Y_{i+1} | y_1^0,...,y_i^0)
    S.C[i+1]<-((new.75[i+1]-C.C[i,(i+1)])/qt(0.75,i+d))^2 
    
    #L_{i+1}
    L.i.plus.1<-matrix(0,i,1) 
    
    for(j in 1:i){
        L.i.plus.1[j,1]<-C.C[j,(i+1)]-y5[(i+1),1]
        }
    
    # Matrix used to calulate T_{i+1}
    M.i.plus.1<-matrix(-y5[1:i],i,i,byrow=T)
    
    if(i>1){
        for(j in 1:(i-1)){
            for(k in (j+1):i){
                M.i.plus.1[j,k]<-M.i.plus.1[j,k] + C.C[j,k]
                }
            }
        }
    for(k in 1:i){
        for(j in k:i){
            M.i.plus.1[j,k]<- M.i.plus.1[j,k]+y0[k]
            }
        }
    
    #T_{i+1}    
    T.i.plus.1<-solve(M.i.plus.1)%*%L.i.plus.1 
    
    # Element (i+1,i+1) of U_{i+1}
    S.i.plus.1<-S.C[i+1]*(1+i/d) / (1 + t(y0[1:i]-y5[1:i])%*%solve(U.i)%*%(y0[1:i]-y5[1:i])/d)  +t(T.i.plus.1)%*%U.i%*%T.i.plus.1
    
    # Construct U_{i+1}
    U.i.plus.1<-matrix(0,i+1,i+1)
    U.i.plus.1[1:i,1:i]<-U.i
    U.i.plus.1[1:i,(i+1)]<-U.i%*%T.i.plus.1
    U.i.plus.1[(i+1),1:i]<-t(U.i.plus.1[1:i,(i+1)])
    U.i.plus.1[(i+1),(i+1)]<-S.i.plus.1
    
    U.i<-U.i.plus.1
    
    # Get estimate of omega
    #F_{i+1}
    F.i.plus.1<-1/(1+ t(y0[1:i]-X[1:i,]%*%b)%*%solve(U.i.plus.1[1:i,1:i])%*%(y0[1:i]-X[1:i,]%*%b)/d)
        
    #H_{i+1}
    H.i.plus.1<-(C.C[(i+1),(i+1)]-C.C[i,(i+1)])/(y0[i+1]-C.C[i,(i+1)])*S.C[i+1]
     
    # Estimate of omega
    omega.hat[i+1]<-(S.C[i+1] - H.i.plus.1 )*(d+i)/d*F.i.plus.1 
    
    }

# Smallest eigenvalue of U_6
min.eig<-min(eigen(U.i)$values)

# Either estimate omega from
#omega<-mean(omega.hat[(omega.hat<min.eig)])

# or specify omega directly
omega<-1085

# Get R.inv
A<-U.i-omega*diag(6)
Xr<-X[1:6,]
R.inv<-solve(t(Xr)%*%Xr)%*%t(Xr)%*%A%*%Xr%*%solve(t(Xr)%*%Xr)/(omega/(d+6))