########################################################################
# Load necessary libraries
########################################################################

library(BMA)
library(MASS)
require(foreign)
require(rrcov)

########################################################################

bicreg.se <- function (x, y, wt = rep(1, length(y)), strict = FALSE, OR = 20, 
    maxCol = 52, drop.factor.levels = TRUE, nbest = 10) 
{
    dropcols <- function(x, y, wt, maxCols = 31) {
        x1.ldf <- data.frame(x, y = y)
        temp.wt <- wt
        lm.out <- lm(y ~ ., data = x1.ldf, weights = temp.wt)
        form.vars <- all.vars(formula(lm.out))[-1]
        any.dropped <- FALSE
        dropped.which <- NULL
        while (length(lm.out$coefficients) > maxCol) {
            any.dropped <- TRUE
            droplm <- drop1(lm.out, test = "none")
            dropped <- row.names(droplm)[which.min(droplm$RSS[-1]) + 
                1]
            dropped.index <- match(dropped, form.vars)
            form.vars <- form.vars[-dropped.index]
            formla <- formula(paste("y", "~", paste(form.vars, 
                collapse = " + "), sep = " "))
            lm.out <- lm(formla, data = x1.ldf, weights = temp.wt)
            dropped.which <- c(dropped.which, dropped)
        }
        new.var.names <- names(lm.out$coefficients)
        return(list(mm = model.matrix(lm.out)[, -1, drop = FALSE], 
            any.dropped = any.dropped, dropped = dropped.which, 
            var.names = new.var.names))
    }
    cl <- match.call()
    x <- data.frame(x)
    if (is.null(dimnames(x))) 
        dimnames(x) <- list(NULL, paste("X", 1:ncol(x), sep = ""))
    y <- as.numeric(y)
    options(contrasts = c("contr.treatment", "contr.treatment"))
    xnames <- dimnames(x)[[2]]
    x2 <- na.omit(data.frame(x))
    used <- match(row.names(data.frame(x)), row.names(x2))
    omitted <- seq(nrow(x))[is.na(used)]
    if (length(omitted) > 0) {
        wt <- wt[-omitted]
        x <- x2
        y <- y[-omitted]
        warning(paste("There were ", length(omitted), "records deleted due to NA's"))
    }
    if (drop.factor.levels) {
        cdf <- cbind.data.frame(y = y, x)
        mm <- model.matrix(formula(cdf), data = cdf)[, -1, drop = FALSE]
        x <- mm
    }
    xx <- dropcols(x, y, wt, maxCol)
    xnames <- xx$var.names[-1]
    x <- xx$mm
    reduced <- xx$any.dropped
    dropped <- NULL
    if (reduced) 
        dropped <- xx$dropped
    nvar <- length(x[1, ])
    if (nvar > 2) {
        a <- leaps(x, y, wt = wt, method = "r2", names = dimnames(x)[[2]], 
            strictly.compatible = FALSE, nbest = nbest)
        a$r2 <- pmin(pmax(0, a$r2), 0.999)
        x.lm <- cbind.data.frame(y = y, as.data.frame(x[, a$which[2, 
            , drop = FALSE]]), w = wt)
        lm.fix <- lm(y ~ . - w, weights = w, data = x.lm)
        r2.fix <- summary(lm.fix)$r.sq
        N <- ncol(x)
        magic <- N * log(1 - a$r2[2]) - N * log(1 - r2.fix)
        a$r2 <- 1 - (1 - a$r2) * exp(-magic/N)
        r2 <- round(c(0, a$r2) * 100, 3)
        size <- c(1, a$size)
        which <- rbind(rep(FALSE, ncol(x)), a$which)
        templabs <- t(matrix(rep(colnames(which), times = nrow(which)), 
            ncol = nrow(which)))
        templabs[!which] <- ""
        label <- apply(templabs, 1, paste, collapse = "")
        label[1] <- "NULL"
    }
    else {
        r2 <- bic <- NULL
        nmod <- switch(ncol(x), 2, 4)
        bic <- label <- rep(0, nmod)
        model.fits <- as.list(rep(0, nmod))
        which <- matrix(c(FALSE, TRUE, FALSE, TRUE, FALSE, FALSE, 
            TRUE, TRUE), nmod, nmod/2)
        size <- c(1, 2, 2, 3)[1:nmod]
        sep <- if (all(nchar(dimnames(x)[[2]]) == 1)) 
            ""
        else ","
        for (k in 1:nmod) {
            if (k == 1) {
                label[k] <- "NULL"
                lm1 <- lm(y ~ 1, w = wt)
            }
            else {
                label[k] <- paste(dimnames(x)[[2]][which[k, ]], 
                  collapse = sep)
                x.lm <- cbind.data.frame(y = y, x = x[, which[k, 
                  , drop = FALSE]], wt = wt)
                lm1 <- lm(y ~ . - wt, data = x.lm, w = wt)
            }
            r2[k] <- summary(lm1)$r.sq * 100
        }
    }
    n <- length(y)
    bic <- n * log(1 - r2/100) + (size - 1) * log(n)
    occam <- bic - min(bic) < 2 * log(OR)
    r2 <- r2[occam]
    size <- size[occam]
    label <- label[occam]
    which <- which[occam, , drop = FALSE]
    bic <- bic[occam]
    mbic <- bic - max(bic)
    postprob <- exp(-0.5 * mbic)/sum(exp(-0.5 * mbic))
    postprob[is.na(postprob)] <- 1
    order.bic <- order(bic, size, label)
    r2 <- r2[order.bic]
    size <- size[order.bic]
    label <- label[order.bic]
    which <- which[order.bic, , drop = FALSE]
    bic <- bic[order.bic]
    postprob <- postprob[order.bic]
    if (strict) {
        nmod <- length(bic)
        if (nmod > 1) {
            occam <- rep(TRUE, nmod)
            for (k in (2:nmod)) {
                for (j in (1:(k - 1))) {
                  which.diff <- which[k, ] - which[j, ]
                  if (all(which.diff >= 0)) 
                    occam[k] <- FALSE
                }
            }
            r2 <- r2[occam]
            size <- size[occam]
            label <- label[occam]
            nmod <- sum(occam)
            which <- which[occam, , drop = FALSE]
            bic <- bic[occam]
            postprob <- postprob[occam]
            postprob <- postprob/sum(postprob)
        }
    }
    probne0 <- round(100 * t(which) %*% as.matrix(postprob), 
        1)
    nmod <- length(bic)
    model.fits <- as.list(rep(0, nmod))
    for (k in (1:nmod)) {
        if (sum(which[k, ]) != 0) {
            model.fits[[k]] <- ls.print(lsfit(x[, which[k, ], 
                drop = FALSE], y, wt = wt), print.it = FALSE)$coef.table[[1]]
        }
        else model.fits[[k]] <- ls.print(lsfit(rep(1, length(y)), 
            y, wt = wt, int = FALSE), print.it = FALSE)$coef.table[[1]]
    }
    Ebi <- rep(0, (nvar + 1))
    SDbi <- rep(0, (nvar + 1))
    CEbi <- Ebi
    CSDbi <- SDbi
    EbiMk <- matrix(rep(0, nmod * (nvar + 1)), nrow = nmod)
    sebiMk <- matrix(rep(0, nmod * (nvar + 1)), nrow = nmod)
    for (i in 1:(nvar + 1)) {
        if ((i == 1) || (sum(which[, (i - 1)] != 0))) {
            for (k in (1:nmod)) {
                if ((i == 1) || (which[k, (i - 1)] == TRUE)) {
                  if (i == 1) 
                    pos <- 1
                  else pos <- 1 + sum(which[k, (1:(i - 1))])
                  EbiMk[k, i] <- model.fits[[k]][pos, 1]

###########################################################
# change x in 
# sebiMk[k, i] <- model.fits[[k]][pos, 2]*x
# to the desired se adjustment (nobs/(nobs-nFEs)) = (37983/(37983-740)) = 1.01986950567892
# below x=2       
###########################################################

                  sebiMk[k, i] <- model.fits[[k]][pos, 2]*1.01986950567892
                }
            }
            Ebi[i] <- as.numeric(sum(postprob * EbiMk[, i]))
            SDbi[i] <- sqrt(postprob %*% (sebiMk[, i]^2) + postprob %*% 
                ((EbiMk[, i] - Ebi[i])^2))
            if (i == 1) {
                CEbi[i] <- Ebi[i]
                CSDbi[i] <- SDbi[i]
            }
            else {
                sel <- which[, i - 1]
                cpp <- postprob[sel]/sum(postprob[sel])
                CEbi[i] <- as.numeric(sum(cpp * EbiMk[sel, i]))
                CSDbi[i] <- sqrt(cpp %*% (sebiMk[sel, i]^2) + 
                  cpp %*% ((EbiMk[sel, i] - CEbi[i])^2))
            }
        }
    }
    dimnames(which) <- list(NULL, dimnames(x)[[2]])
    dimnames(EbiMk) <- dimnames(sebiMk) <- list(NULL, c("Int", 
        dimnames(x)[[2]]))
    result <- list(postprob = postprob, namesx = xnames, label = label, 
        r2 = r2, bic = bic, size = (size - 1), which = which, 
        probne0 = c(probne0), postmean = Ebi, postsd = SDbi, 
        condpostmean = CEbi, condpostsd = CSDbi, ols = EbiMk, 
        mle = EbiMk, se = sebiMk, reduced = reduced, dropped = dropped, 
        call = cl, n.models = length(postprob), n.vars = length(probne0))
    class(result) <- "bicreg"
    result
}


###########################################################################################
#iterative
###########################################################################################
iterative.bicreg.se <- function (x, y, maxVars,
        threshold = 1, adaptive.threshold =1, wt = rep(1, length(y)), strict = FALSE, OR = 20,
       maxCol = maxVars+1, drop.factor.levels = TRUE, nbest = 10){

   ## consider maxVars at a time in traditional bma
   ## drop variables with probne0 < threhsold / probne0 <  (min(probne0) +threshold)

   ### rank X matrix
   temp <- NULL
   for(i in 1:ncol(x))
    temp <- c(temp , summary(lm(y~x[,i]))$r.squared)
   x2 <- x[, rev(order(temp))]

   this.x <- x2[ , 1:maxVars ]
   x2 <- x2[ , -c(1:maxVars )]

   loop<-function(){
     bicreg.call <-function() bicreg.se(this.x, y, wt = wt, strict = strict, OR=OR, maxCol = maxCol, drop.factor.levels=drop.factor.levels, nbest=nbest)

     # while we have variables left to be processed.
     adapt.dropped <- NULL
     while(ncol(x2)>0){

     # traditional bma
     model <- bicreg.call()

     # drop varaibles
     dropme <- which(model$probne0 < threshold  )
     # adaptive threshold
     if(length(dropme)==0){
         dropme <- which(model$probne0 < min(model$probne0 + adaptive.threshold))
         print(paste(model$namesx[dropme], "with probne0 = ", model$probne0[dropme], "dropped in adaptive threshold step"))
	   adapt.dropped <- c(adapt.dropped, model$namesx[dropme])
     }

     # update variables to consider
     this.x<-this.x[ , -c(dropme)]
     if(ncol(x2)<length(dropme)) dropme<-1:ncol(x2)

     old.names.x2 <- names(x2)
     newdata<-data.frame(x2[, 1:length(dropme)])
     names(newdata)<-old.names.x2[1:length(dropme)]
     this.x<-cbind(this.x, newdata)
     x2<-data.frame(x2[, -c(1:length(dropme))])
     names(x2)<-old.names.x2[-c(1:length(dropme))]

  }
  assign("adapt.dropped", adapt.dropped, .GlobalEnv)
  bicreg.call()
}

  # traditional bma, then reverse loop on vars dropped in adaptive threshold
  this.x<-x[,loop()$namesx]
  x2<-x[,rev(adapt.dropped)]
  print("---------------------------------------------")
  loop()
}



#****************************************************************
#Basic Variables for tvc-Dataset (23)
#****************************************************************

# Trade Creation: (13)
# MX_asean92, MX_anzcerta, MX_apec, MX_ap, MX_cacm, MX_caricom, MX_eea, MX_efta, MX_eu, MX_laia, MX_mercosur, MX_nafta, MX_bi,

# Core Gravity: (1)
# ldist,

# Economic Policy: (2)
# custrict, vola3,

# Development, Factor Endowment: (3)
# diffgdppc, diffdensity, diffsecschool25,

# Geography: (1)
# border,

# History: (3)
# comlang, comcol, colony


#****************************************************************
#Advanced Variables for tvc-Dataset
#****************************************************************
# OVERLAPPING
 
# Overlapping Trade Creation - ONE period BEFORE accession: (13)
# MXm1_asean92, MXm1_anzcerta, MXm1_apec, MXm1_ap, MXm1_cacm, MXm1_caricom, MXm1_eea, MXm1_efta, MXm1_eu, MXm1_laia, MXm1_mercosur, MXm1_nafta, MXm1_bi,

# Overlapping Trade Creation - ONE period AFTER accession: (13)
# MXp1_asean92, MXp1_anzcerta, MXp1_apec, MXp1_ap, MXp1_cacm, MXp1_caricom, MXp1_eea, MXp1_efta, MXp1_eu, MXp1_laia, MXp1_mercosur, MXp1_nafta, MXp1_bi,

# Overlapping Trade Creation - TWO period AFTER accession: (10)
# MXp2_anzcerta, MXp2_apec, MXp2_ap, MXp2_cacm, MXp2_caricom, MXp2_efta, MXp2_eu, MXp2_laia, MXp2_nafta, MXp2_bi,

# Overlapping Trade Creation - one period AFTER DROPPING OUT: (2)
# MXa1_efta, MXa1_ap


#****************************************************************
# NON-OVERLAPPING

# Non-overlapping Trade Creation - at time of accession: (13)
# MXo_asean92, MXo_anzcerta, MXo_apec, MXo_ap, MXo_cacm, MXo_caricom, MXo_eea, MXo_efta, MXo_eu, MXo_laia, MXo_mercosur, MXo_nafta, MXo_bi,

# Non-overlapping Trade Creation - ONE period BEFORE accession: (13)
# MXom1_asean92, MXom1_anzcerta, MXom1_apec, MXom1_ap, MXom1_cacm, MXom1_caricom, MXom1_eea, MXom1_efta, MXom1_eu, MXom1_laia, MXom1_mercosur, MXom1_nafta, MXom1_bi,

# Non-overlapping Trade Creation - ONE period AFTER accession: (10)
# MXop1_anzcerta, MXop1_apec, MXop1_ap, MXop1_cacm, MXop1_caricom, MXop1_efta, MXop1_eu, MXop1_laia, MXop1_nafta, MXop1_bi,

# Non-overlapping Trade Creation - one period AFTER DROPPING OUT: (2)
# MXoa1_efta, MXoa1_ap


#****************************************************************
# CONTROL VARIABLES BROKEN DOWN BY EX/IM:

# Core Gravity: (1)
# ldist, 

# Economic Policy: (2)
# custrict, vola3,

# Development, Factor Endowment: (3)
# diffgdppc, diffdensity, diffsecschool25,

# Geography: (1)
# border,

# History: (3)
# comlang, comcol, colony


#****************************************************************
# OTHER VARIABLES:

# Variables multiplied together, but take square root:
# lgdpij, lgdppcij, lremoteij, lareaij,

# Real GDP:
# lrgdp1, lrgdp2, lrgdpij, lprdctrgdpij,

# Other variables not in GY:
# curcol, M_WTO, X_WTO, MX_WTO, TD_WTO


#############################################################################
# tvc Regressions
#############################################################################

setwd("C:\\Documents and Settings\\te\\Desktop\\christian")
load("tvc.RData")


y <- labsimport

x3 <- cbind(MXom1_asean92, MXom1_anzcerta, MXom1_apec, MXom1_ap, MXom1_cacm, MXom1_caricom, MXom1_eea, MXom1_efta, MXom1_eu, MXom1_laia, MXom1_mercosur, MXom1_nafta, MXom1_bi, MXo_asean92, MXo_anzcerta, MXo_apec, MXo_ap, MXo_cacm, MXo_caricom, MXo_eea, MXo_efta, MXo_eu, MXo_laia, MXo_mercosur, MXo_nafta, MXo_bi, MXp1_asean92, MXp1_anzcerta, MXp1_apec, MXp1_ap, MXp1_cacm, MXp1_caricom, MXp1_eea, MXp1_efta, MXp1_eu, MXp1_laia, MXp1_mercosur, MXp1_nafta, MXp1_bi, MXoa1_efta, MXoa1_ap, ldist, custrict, vola3, diffgdppc, diffdensity, diffsecschool25, border, comlang, comcol, colony)

x3 <- as.matrix(x3)

# Iterative Bicreg with DF adjustment fcn is:
# iterative.bicreg.se <- function (x, y, maxVars, threshold = 0.01, adaptive.threshold =1, wt = rep(1, length(y)), strict = FALSE, OR = 20, maxCol = maxVars+1, drop.factor.levels = TRUE, nbest = 10)

ibma_tvc_res3 <- iterative.bicreg.se(x=x3, y=y, maxVars=36 , OR = 20, drop.factor.levels = TRUE, nbest=10)
# OR=10^20 nbest=30
save.image("c:\\aaa_BMA\\tvc.RData")

summary(ibma_tvc_res3, conditional=TRUE)


save.image("c:\\aaa_BMA\\tvc.RData")



# no breakdown by im-ex is needed for tvc, as all im & ex-specific variables are dropped anyways