Holt Winters Forecast with Multiple Input Variables

Question

For context, I'm a novice R user, so please forgive any incorrect terminology/processes. I am actively trying to improve my coding ability, but recently have become stumped.

I have the following data set where A * B * C = Output:

Date          A           B                     C                     Output
1/1/2013   177352 0.908329198 0.237047935 38187
1/2/2013   240724 0.852033865 0.237273592 48666
1/3/2013   243932 0.908380204 0.237039845 52524
1/4/2013   221485 0.820543152 0.236356733 42955
1/5/2013   202590 0.818066045 0.240900973 39925
1/6/2013   238038 0.770057722 0.247344561 45339
1/7/2013   271511 0.794258796 0.241252029 52026
1/8/2013   283434 0.807817693 0.233810703 53534
1/9/2013   275016 0.843220031 0.243769917 56530
1/10/2013 255266 0.797791324 0.238562428 48583
1/11/2013 226564 0.815791564 0.236153417 43648
1/12/2013 214366 0.800066242 0.237961133 40812
1/13/2013 256946 0.764845532 0.237640186 46702
1/14/2013 282298 0.816537843 0.234257528 53998

I have a few years worth of data and I'm trying for forecast Output, using A, B, and C. However, when I model out A, B, and C individually, the Output becomes very skewed. If I forecast just Output then I lose the input factors.

What is the best package/code to accomplish this task? I've tried Googling and searching on here numerous different methods, but haven't found the solution I'm looking for.

Here is some of the code:

DataSet1[,"Date"] <- mdy(DataSet[,"Date"])

DataSet1
TotalSet <- ts(DataSet1, frequency = 365, start =c(2013,1))
DataA <- ts(DataSet1$A, frequency = 365, start = c(2013,1))
DataB <- ts(DataSet1$B, frequency = 365, start = c(2013,1))
DataC <- ts(DataSet1$C, frequency = 365, start = c(2013,1))
OutputData <- ts(DataSet$Output, frequency = 365, start = c(2013,1))

ADecompose <- decompose(DataA)
BDecompose <- decompose(DataB)
CDecompose <- decompose(DataC)
OutputDecompose <- decompose(OutputData)

DataAHW <- HoltWinters(DataA, seasonal = "mult")
DataBHW <- HoltWinters(DataB, seasonal = "mult")
DataCHW <- HoltWinters(DataC, seasonal = "mult")
OutputDataHW <- HoltWinters(OutputData, seasonal = "mult")

FC.A <- forecast.HoltWinters(DataAHW)
FC.B <- forecast.HoltWinters(DataBHW)
FC.C <- forecast.HoltWinters(DataCHW)
FC.Output <- forecast.HoltWinters(OutputDataHW)

plot(ForecastVisits)
plot(ForecastCPV)
plot(ForecastRPC)
plot(ForecastRevenue)

Here is another model I built for the Output and I've plugged A, B, and C into it individually then combined them in excel. I'm sure there is a more appropriate way to handle this, but given my lack of experience I am reaching out for help

dataset <- testData


##FORECAST

forecastingFuntion <- function(dataset, lenghtOfForecast)

  {
dataset[,"Date"] <- mdy(dataset[,"Date"])
myts <- ts(dataset[,"DataSet$Output"], start = c(2013,1), frequency = 365)

hwModel <- HoltWinters(myts, seasonal = "mult")

future <- data.frame(predict(hwModel, n.ahead = lenghtOfForecast, level =       0.9))
fittedValues <- data.frame(as.numeric(hwModel$fitted[,"xhat"]))
names(fittedValues) <- "fit"

futureDates <- c()
predicitedValues <- rbind(fittedValues, future)

for(i in 1: lenghtOfForecast)
  {
  futureDateSingle <- data.frame(dataset[nrow(dataset),"Date"] + days(i))
  futureDates <- rbind(futureDates, futureDateSingle)  
  }  

names(futureDates) <- "Date"
dates <- data.frame(dataset[366:(nrow(dataset)),"Date"])
names(dates) <- "Date"
dates <- rbind(dates, futureDates)
predictedData <- data.frame(predicitedValues, dates)
names(predictedData) <- c("predictedValues","Date")
finalData2 <- mergeData <- merge(predictedData, dataset, all.x = T, all.y = F, by = "Date")
finalData2
}



finalData2 <- forecastingFuntion(testData, 612)

rm(list=setdiff(ls(), c("finalData2")))

write.csv(finalData2, file="B2BForecastVisits.csv")

Thanks!