SalaiMartinFrischWaugh <- function(dependent, regressors, examined_regressor,keepfix_vector,include_intercept=TRUE,number_of_changing_regressors=3) {
    
    
    yall=data.matrix(dependent)
    Xall=data.matrix(regressors)
    yXall=na.omit(cbind(yall,Xall))
    ncolall=ncol(yXall)-1
    
    #usability checks
    examined_regressor=examined_regressor[1]
    if (is.character(examined_regressor)) {
        examined_regressor=which(names(regressors)==examined_regressor)
    }
    if (is.character(keepfix_vector)) {
        keepfix_vector=which(!(is.na(match(names(regressors),keepfix_vector))))
    }

    examined_regressor=as.integer(examined_regressor)
    keepfix_vector=as.integer(keepfix_vector)


    
    #do frisch-waugh for variables kept fixed
    if (sum(keepfix_vector)>0) {
        #if there are variables to be kept fixed, 
        #condition the others on them by applying Frisch-Waugh
        FWdesign=yXall[,(keepfix_vector+1)]
        if (include_intercept) {
            FWdesign=cbind(1,FWdesign)
        }
        FWregressors=yXall[,-(keepfix_vector+1)]
        FWbeta=chol2inv(chol(crossprod(FWdesign))) %*% crossprod(FWdesign,FWregressors)
        yXall_cleaned=FWregressors-FWdesign %*% FWbeta
        ncolall_cleaned=ncolall-length(keepfix_vector)
        examined_index=which((1:ncolall)[-keepfix_vector]==examined_regressor)

    } else {
        #if there are no variables to be kept fixed
        ncolall_cleaned=ncolall
        examined_index=examined_regressor
        if (include_intercept) {
            #using Frisch-Waugh for the incerpt: demean everything
            yXall_cleaned = yXall - apply(yXall,MARGIN=2,FUN="mean")
        } else {
            yXall_cleaned=yXall
        }

    }
        

    #get combinations
    combiset=combn((1:ncolall_cleaned)[-examined_index],number_of_changing_regressors)
     
    #initialize storing variables
        total_lik <- 0
        total_beta <- 0
        total_betavariance <- 0
        total_betavarianceWhite <- 0
        totalW_beta <- 0
        totalW_betavariance <- 0
        totalW_betavarianceWhite <- 0


    cat(as.character(Sys.time()),": Starting Simulation runs...\n")

    for (irun in 1:ncol(combiset)) {
        regindices2include=combiset[,irun]
        yX2regress=yXall_cleaned[,c(1,examined_index+1,(regindices2include+1))]
    
        model=OLSest(yXmat=yX2regress,bincludeintercept=FALSE,bwhite=TRUE)
    
        #sum up 'on the fly' and weight with likelihood
        currlik=exp(model$loglik)
        totalW_beta<-totalW_beta + model$betahat[1]*currlik
        totalW_betavariance <- totalW_betavariance + as.vector(model$betavar[1,1])*currlik
        totalW_betavarianceWhite <- totalW_betavarianceWhite + as.vector(model$betavar_white[1,1])*currlik
        total_lik <- total_lik + currlik

        total_beta<-total_beta + model$betahat[1]
        total_betavariance <- total_betavariance + as.vector(model$betavar[1,1])
        total_betavarianceWhite <- total_betavarianceWhite + as.vector(model$betavar_white[1,1])

    
        if (irun %% 5000==0) {cat(as.character(Sys.time()),": run number ",irun,"\n")}
    }
    
    av_beta <- total_beta / irun
    av_betavariance <- total_betavariance / irun
    av_betaSE <- av_betavariance^0.5
    av_betavarianceWhite <- total_betavarianceWhite / irun
    av_betaWhiteSE <- av_betavarianceWhite^0.5
    pval_t <- (1-pt(abs(av_beta/av_betaSE),nrow(yXall_cleaned)))*2
    pval_t_white <- (1-pt(abs(av_beta/av_betaWhiteSE),nrow(yXall_cleaned)))*2
    
    Wav_beta <- totalW_beta / total_lik
    Wav_betavariance <- totalW_betavariance / total_lik
    Wav_betaSE <- Wav_betavariance^0.5
    Wav_betavarianceWhite <- totalW_betavarianceWhite / total_lik
    Wav_betaWhiteSE <- Wav_betavarianceWhite^0.5
    Wpval_t <- (1-pt(abs(Wav_beta/Wav_betaSE),nrow(yXall_cleaned)))*2
    Wpval_t_white <- (1-pt(abs(Wav_beta/Wav_betaWhiteSE),nrow(yXall_cleaned)))*2

    #Output
    
    av_outlist=list(betahat=av_beta ,betaSE=av_betaSE,Pval=pval_t,betaSEwhite=av_betaWhiteSE,Pval_white=pval_t_white)
    Wav_outlist=list(betahat=Wav_beta ,betaSE=Wav_betaSE,Pval=Wpval_t,betaSEwhite=Wav_betaWhiteSE,Pval_white=Wpval_t_white)

    #Display nice information
    av_outtable=av_outlist
    names(av_outtable)=c("Av.Beta","Std.Err.","Prob.","White Std.Err.","White Prob.")
    Wav_outtable=Wav_outlist
    names(Wav_outtable)=c("Av.Beta","Std.Err.","Prob.","White Std.Err.","White Prob.")

    cat("\n\n--------- RESULTS ------------------------------------\n")
    cat("   Tested variable number",examined_regressor,":",names(regressors)[examined_regressor],"\n")
    if (sum(keepfix_vector)>0) {
        cat("   Fixed regressors",paste(keepfix_vector, collapse=", "),":",names(regressors)[keepfix_vector],"\n")
    } else {
        cat("   Fixed regressors: none\n")
    }
    cat("  ",irun,"iterations over all other regressors\n")
    cat("\n------- WEIGHTED AVERAGE PARAMETERS -------------------------------------------\n")
    print(unlist(Wav_outtable))
    cat("\n------- SIMPLE AVERAGE PARAMETERS ----------------------------------------------\n")
    print(unlist(av_outtable))
    cat("--------------------------------------------------------------------------------\n")

    return(list(weighted=Wav_outlist, simple=av_outlist))
}

#millions.csv ist just the excel file millions.xls saved as CSV (in comma-separated-value format), nothing else

growthdata = read.csv("millions.csv",na.strings=c("#NV","."),header=T,row.names=2,nrows=134)

#we have to tell R to read in only 134 rows, since the excel file contained several more rows of junk
#moreover we set row.names=2 (the country codes are in the second column), and specify what NA is in the Excel CSV format.



testoutput1=SalaiMartinFrischWaugh(growthdata[,3],growthdata[,4:64],"X39",c("X1","X2","X3"),T)