Ok, so here is what I have to offer that may be of use to some people in my quest for donations:
There are two R scripts, the first one downloads all the trades from bitcoincharts (trades are aggregated together if they occur within a second of each other I believe) and saves it as a csv. After this is done the first time, downloading all the trades will be skipped and it will only download trades since the last time the script was run and update the csv file. It then calculates the daily data from that.
It next gets the bitcointalk forum stats as well as wikipedia page views and makes a chart shown below. Then what it does is take the price history, calculate "breakpoints" which indicate timepoints of structural change and overlays a plot of this onto a heatmap showing possible lines of resistance according to how much volume has occurred at each price. Volume can be in USD, BTC, or # of trades per day. This is the second plot (using USD volume).
The second script will take this daily data, and attempt to predict into the future based on previous day-to-day price changes as has been shown throughout this thread. It can also take snapshots of what it is doing and save them to a folder which can then be made into a video. I recommend using virtual dub for this.
It was made on windows 7 using R 2.14.2 (64 bit), but should work using the newer versions of R. Not sure what problems will arise on other OSes.
Script 1:
####Only Run First Time (uncomment to run)
#install.packages("rjson")
#install.packages("TTR")
#install.packages("RCurl")
#install.packages("XML")
#install.packages("fields")
#install.packages("psych")
#install.packages("strucchange")
#Load required packages
require(rjson)
require(TTR)
require(RCurl)
require(XML)
require(fields)
require(psych)
require(strucchange)
#Choose Directory for data
data.dir<-"C:/Users/Valued Customer/Documents/GoxData"
# Download All gox USD history from bitcoincharts.com? T or F
# If false it will only update from previous file in data directory named "goxUSD.csv"
dl.all=F
# Get Google Trends Data (experimental, doesnt work)
get.gtrends=F
#Plot Breakpoints and volume heatmaps? This is slow.
plot.breakpoints=T
#Volume for heatmaps: "usd", "trades", or "btc"
heatmap.vol="usd"
#Set working Directory
setwd(data.dir)
####Experimental google trends (requires copy/paste from exported csv)####
if(get.gtrends==T){
new.gtrends=F
if(new.gtrends==T){
gtrends1<-readClipboard()
out=NULL
for(i in 1:length(gtrends1)){
temp<-strsplit(gtrends1[i]," ")[[1]][1]
out<-rbind(out,temp)
}
gtrends1<-out
out=NULL
for(i in 1:nrow(gtrends1)){
date<-strsplit(as.character(gtrends1[i,1]), "-")
yr<-as.numeric(date[[1]][1])
mo<-as.numeric(date[[1]][2])
day<-as.numeric(date[[1]][3])
temp<-cbind(yr,mo,day)
out<-rbind(out,temp)
}
gtrends<-data.frame()
gtrends<-as.data.frame(cbind(out,c(gtrends1),gtrends2))
gtrends[,1]<-as.numeric(as.character(gtrends[,1]))
gtrends[,2]<-as.numeric(as.character(gtrends[,2]))
gtrends[,3]<-as.numeric(as.character(gtrends[,3]))
gtrends[,5]<-as.numeric(as.character(gtrends[,5]))
gtrends<-gtrends[which(gtrends[,4]=="2009-11-29"):nrow(gtrends),]
gtrends<-gtrends[is.finite(gtrends[,5]),]
colnames(gtrends)<-c("Year","Month","Day","Date","Value")
rownames(gtrends)<-1:nrow(gtrends)
write.table(gtrends,"gtrends.txt")
}else{
gtrends<-read.table("gtrends.txt")
}
}
##End Google Trends
####Get All goxUSD data###
if( dl.all==T){
gox<-read.table("http://api.bitcoincharts.com/v1/csv/mtgoxUSD.csv",sep=",")
##Convert Times from Unix
times=matrix(nrow=nrow(gox), ncol=6)
time.string=matrix(nrow=nrow(gox),ncol=1)
pb<-txtProgressBar(min = 0, max = nrow(gox), initial = 0, char = "=",
width = NA, style = 3)
for(i in 1:nrow(gox)){
timeStamp<-ISOdatetime(1970,1,1,0,0,0) + gox[i,1]
date<-strsplit(as.character(timeStamp), " ")[[1]][1]
yr<-as.numeric(strsplit(date, "-")[[1]][1])
mo<-as.numeric(strsplit(date, "-")[[1]][2])
day<-as.numeric(strsplit(date, "-")[[1]][3])
ToD<-strsplit(as.character(timeStamp), " ")[[1]][2]
hr<-as.numeric(strsplit(ToD, ":")[[1]][1])
min<-as.numeric(strsplit(ToD, ":")[[1]][2])
sec<-as.numeric(strsplit(ToD, ":")[[1]][3])
times[i,]<-cbind(yr,mo,day,hr,min,sec)
time.string[i,]<-as.character(timeStamp)
setTxtProgressBar(pb, i)
}
close(pb)
gox<-cbind(times,gox)
gox<-cbind(gox,gox[,8]*gox[,9])
colnames(gox)<-c("Yr","Mo","Day","Hr","Min","Sec","UnixT","Price","BTCVol", "USDVol")
write.csv(gox,"goxUSD.csv", row.names=F)
}
##End get All goxUSD
####Update goxUSD data with new prices####
gox<-read.csv("goxUSD.csv", header=T, sep=",")
while((tail(gox[,7],1)+60*15)
update<-read.table(paste("http://api.bitcoincharts.com/v1/trades.csv?symbol=mtgoxUSD&start=",
(tail(gox[,7],1)+1),sep="")
,sep=",")
##Convert Times from Unix
times=matrix(nrow=nrow(update), ncol=6)
time.string=matrix(nrow=nrow(update),ncol=1)
pb<-txtProgressBar(min = 0, max = nrow(update), initial = 0, char = "=",
width = NA, style = 3)
for(i in 1:nrow(update)){
timeStamp<-ISOdatetime(1970,1,1,0,0,0) + update[i,1]
date<-strsplit(as.character(timeStamp), " ")[[1]][1]
yr<-as.numeric(strsplit(date, "-")[[1]][1])
mo<-as.numeric(strsplit(date, "-")[[1]][2])
day<-as.numeric(strsplit(date, "-")[[1]][3])
ToD<-strsplit(as.character(timeStamp), " ")[[1]][2]
hr<-as.numeric(strsplit(ToD, ":")[[1]][1])
min<-as.numeric(strsplit(ToD, ":")[[1]][2])
sec<-as.numeric(strsplit(ToD, ":")[[1]][3])
times[i,]<-cbind(yr,mo,day,hr,min,sec)
time.string[i,]<-as.character(timeStamp)
setTxtProgressBar(pb, i)
}
close(pb)
update<-cbind(times,update)
update<-cbind(update,update[,8]*update[,9])
colnames(update)<-c("Yr","Mo","Day","Hr","Min","Sec","UnixT","Price","BTCVol", "USDVol")
gox<-rbind(gox,update)
}
#Remove Duplicate Entries (maybe bitcoincharts bug?)
if(length(which(duplicated(gox)))>0){
gox<-gox[-which(duplicated(gox)),]
}
write.csv(gox,"goxUSD.csv", row.names=F)
##End Update goxUSD
####Aggregate Gox data into Daily Info####
Unix.start<-gox[1,7] - gox[1,4]*3600 - gox[1,5]*60 - gox[1,6]
daily=NULL
while(Unix.start
Unix.end<-Unix.start+60*60*24
temp<-gox[which(gox[,7]>Unix.start & gox[,7]
temp2<-cbind(
temp[1,1],
temp[1,2],
temp[1,3],
Unix.start,
sum(temp[,10])/sum(temp[,9]),
sum(temp[,9]),
sum(temp[,10]),
nrow(temp)
)
colnames(temp2)<-c("Yr","Mo","Day","UnixT","VWAP","BTCVol","USDVol","num.Trades")
daily<-rbind(daily,temp2)
Unix.start<-Unix.end + 1
}
first.gox<-paste(daily[1,1],"-",0,daily[1,2],"-",daily[1,3],sep="")
daily[(which(is.na(daily[,5]))),5]<-
daily[(which(is.na(daily[,5]))[1]-1),5]
write.csv(daily,"goxUSDdaily.csv", row.names=F)
##End Daily
##Get months from unix timestamps
month.labels=data.frame()
for(y in 2010:2013){
for(m in 1:12){
unixT<-gox[which(gox[,1]==y & gox[,2]==m),7][1]
month.labels<-rbind(month.labels,cbind(unixT,month.abb[m]))
}
}
month.labels<-month.labels[which(!is.na(month.labels[,1])),]
month.labels<-month.labels[2:nrow(month.labels),]
month.labels[,2]<-as.character(month.labels[,2])
yrs<-c(11,12,13)
for(i in 1:length(which(month.labels[,2]=="Jan"))){
month.labels[which(month.labels[,2]=="Jan"),2][1]<-paste(
month.labels[which(month.labels[,2]=="Jan"),2][1],
yrs[i], sep="")
}
##End Month Labels
####Get Forum Data####
forum2=NULL
for(y in 2009:2013){
for (m in 1:12){
start<-paste(month.name[m],y)
if(m==1){
end<-paste(month.name[12],(y-1))
}else{
end<-paste(month.name[m-1],y)
}
date<-as.numeric(paste(y,m,sep=""))
forum<-getURL(paste("https://bitcointalk.org/index.php?action=stats;expand=",
date,"#",date,sep=""),
ssl.verifypeer = FALSE, useragent = "R"
)
forum<-readLines(tc <- textConnection(forum)); close(tc)
pagetree <- htmlTreeParse(forum, error=function(...){}, useInternalNodes = TRUE)
x <- xpathSApply(pagetree, "//*/table", xmlValue)
x <- unlist(strsplit(x, "\n"))
x <- gsub("\t","",x)
x <- sub("^[[:space:]]*(.*?)[[:space:]]*$", "\\1", x, perl=TRUE)
x <- x[!(x %in% c("", "|"))]
if(any(x==start) & any(x==end)){
temp<-as.data.frame(rbind(x[which(x==start)[2]:(which(x==end)[2]-1)]))
out=NULL
for(j in 2:(length(temp)/6)){
n<-seq(1,length(temp),by=6)[j]
temp2<-as.matrix(temp[n:(n+5)])
out<-rbind(out,temp2)
}
}else{
out<-matrix(rep(NA,6),nrow=1,ncol=6)
}
forum2<-rbind(forum2,out)
}
}
forum2<-forum2[!is.na(forum2[,1]),]
colnames(forum2)<-c("Date","New Topics","New Posts","New Members","Most Online","Page Views")
forum<-as.data.frame(forum2)
for(i in 2:6){
forum[,i]<-as.numeric(as.character(forum[,i]))
}
##End Forum Stats
####Get Wikipedia Page Views####
wiki3=NULL
for(y in 2009:2013){
for (m in 1:12){
if(m<10){
mo<-paste(0,m,sep="")
}else{
mo<-m
}
date<-as.numeric(paste(y,mo,sep=""))
wiki<-fromJSON(
getURL(
paste("http://stats.grok.se/json/en/",
date,
"/Bitcoin",sep="")
)
)
wiki<-as.matrix(wiki$daily_views)
wiki2<-cbind(rownames(wiki),wiki)
rownames(wiki2)<-1:nrow(wiki2)
out=NULL
for(i in 1:nrow(wiki2)){
date2<-strsplit(as.character(wiki2[i,1]), "-")
yr<-as.numeric(date2[[1]][1])
mo<-as.numeric(date2[[1]][2])
day<-as.numeric(date2[[1]][3])
temp<-cbind(yr,mo,day)
out<-rbind(out,temp)
}
wiki2<-cbind(out[,1],out[,2],out[,3],wiki2)
wiki2<-as.data.frame(wiki2)
wiki2[,1]<-as.numeric(wiki2[,1])
wiki2[,2]<-as.numeric(wiki2[,2])
wiki2[,3]<-as.numeric(wiki2[,3])
wiki2[,5]<-as.numeric(wiki2[,5])
wiki2<-wiki2[order(wiki2[,3]),]
rownames(wiki2)<-1:nrow(wiki2)
wiki3<-rbind(wiki3,wiki2)
}
}
shrt.months<-c(2,4,6,9,11)
lp.yrs<-seq(2000,2020,by=4)
del.points<-c(which(wiki3[,2] %in% shrt.months & wiki3[,3]==31),
which(wiki3[,2]==2 & wiki3[,3]==30),
which(!wiki3[,1] %in% lp.yrs &
wiki3[,2]==2 & wiki3[,3]==29)
)
wiki3<-wiki3[-del.points,]
rownames(wiki3)<-1:nrow(wiki3)
wiki<-wiki3[which(wiki3[,4]==forum[1,1]):nrow(wiki3),]
rownames(wiki)<-1:nrow(wiki)
##End Wikipedia Page Views
####Plot Interest Indicators####
x.axis.labels=NULL
for(i in seq(1,nrow(forum),by=7*8)){
x.axis.labels<-rbind(x.axis.labels,paste(forum[i,1]))
}
dev.new()
plot(
forum[,2]/max(forum[,2]),
ylab="Value Normalized to Max",
xaxt="n", xlab="Time",
col=rainbow(6, alpha=.75)[1], type="l", lwd=2,
main="Normalized to Series Max"
)
axis(1,at=seq(1,nrow(wiki),by=7*8),labels=x.axis.labels, las=1, cex=.7)
for(i in 3:5){
lines(
forum[,i]/max(forum[,i]), col=rainbow(6, alpha=.75)[i-1], lwd=2,
)
}
lines(wiki[,5]/max(wiki[,5]), col=rgb(0,0,0,.75), lwd=2)
if(get.gtrends==T){
lines(seq(1,nrow(wiki),by=7),gtrends[1:(nrow(gtrends)),5]/
max(gtrends[1:(nrow(gtrends)-1),5]), col=rainbow(6, alpha=.75)[5],lwd=4
)
}
lines(which(forum[,1]==first.gox):nrow(forum),
daily[,5]/max(daily[,5],na.rm=T),
lwd=6,lty=1,col="Black"
)
lines(which(forum[,1]==first.gox):nrow(forum),
daily[,5]/max(daily[,5],na.rm=T),
lwd=3,lty=2,col="Yellow"
)
legend("topleft", c(paste("BitcoinTalk:",colnames(forum)[2:5]),"Google Trends (Bitcoin + Bitcoins)", "Wikipedia Page Hits", "MtGox USD/BTC Daily VWAP"),
col=c(rainbow(6, alpha=.75),rgb(0,0,0,.75),"Black"), lty=c(rep(1,6),1),
lwd=c(4,4,4,4,6,4,8),bty="n",y.intersp=.75
)
legend("topleft", c(paste("BitcoinTalk:",colnames(forum)[2:5]),"Google Trends (Bitcoin + Bitcoins)", "Wikipedia Page Hits", "MtGox USD/BTC Daily VWAP"),
col=c(rainbow(6, alpha=.75)[1:5],rgb(0,0,0,.75),"Yellow"), lty=c(rep(1,6),2),
lwd=c(4,4,4,4,6,4,5),bty="n",y.intersp=.75
)
##End Plots
if(plot.breakpoints==T){
####Plot Breakpoints with Volume Heatmaps (slow)####
trades=F
usd=T
btc=F
t<-(1:length(daily[,5]))
bp<-breakpoints(daily[,5]~t, h=30)
##Format data##
price <- round(daily[,5],2) # for line plot
if(heatmap.vol=="trades"){
volume <- daily[,8] # for line color heatmap
legend.lab="Trades per Day (Normalized to Max)"
}
if(heatmap.vol=="usd"){
volume <- daily[,7] # for line color heatmap
legend.lab="USD per Day (Normalized to Max)"
}
if(heatmap.vol=="btc"){
volume <- daily[,6] # for line color heatmap
legend.lab="BTC per Day (Normalized to Max)"
}
time <- 1:length(price) # X axis
##Smoothing Kernel Settings##
FWHM<-1 #Set size of smoothing kernal (in pixels)
time.resolution<-.00001
price.resolution<-.05
lower.cutoff.multiplier<-.01
##Create Heatmap Matrix##
price.range<-round(seq(min(price,na.rm=T),(max(price,na.rm=T)+1),by=.1),2)
heatmap<-matrix(nrow=length(time),ncol=length(price.range))
for (i in 1:length(time)){
heatmap[i,which(price.range==price[i])]<-volume[i]
}
heatmap[which(is.na(heatmap))]<-0
##Smooth the heatmap##
smooth.heatmap<-image.smooth(heatmap, dx=time.resolution, dy=price.resolution, Nwidth=, Mwidth=,
theta=FWHM, kernel.function=dnorm, tol= 1.1, na.rm=TRUE)
smooth.heatmap$z[which(smooth.heatmap$z<0)]<-0
smooth.heatmap$z<-smooth.heatmap$z/max(smooth.heatmap$z)
##Heatmap Chart Settings###
max.smoothedvolume<-max(smooth.heatmap$z)
lower.cutoff<-lower.cutoff.multiplier*max.smoothedvolume
number.of.color.bins<-1000
colors<-c("Grey",rev(rainbow(number.of.color.bins-1)))
##Make Chart##
dev.new()
image.plot(time, price.range, as.matrix(smooth.heatmap$z),
breaks= c(0,seq(lower.cutoff, max.smoothedvolume, length=number.of.color.bins)),
col=colors, legend.lab=legend.lab,
ylab="USD/BTC (Mt.Gox VWAP)", xlab="Time"
)
#lines(time, price, lwd=4)
lines(daily[,5], type="l", lwd=3)
points(as.numeric(names(bp[[9]])[bp[[1]]]),
bp[[9]][bp[[1]]],
col="Black",pch=16, type="b", lwd=3,cex=1.2
)
lines(daily[,5], type="l", lwd=3)
points(as.numeric(names(bp[[9]])[bp[[1]]]),
bp[[9]][bp[[1]]],
col="Red",pch=16, type="b"
)
##End Breakpoints and Heatmaps
}
Script 2
require("fields")
#########Settings##############
# Set data directory and directory for movie frames
data.dir<-"C:/Users/Valued Customer/Documents/GoxData"
movie.dir<-"C:/Users/Valued Customer/Documents/GoxPredict/"
days.before.today.start<-0 #When To start?
iters<-1000 #Number of Simulations to Run
pred.days<-30 #Days Ahead to Predict
num.examps<-10 #Number of Example Scenarios to Plot
live.plots=T #Make Movie?
make.pngs=F
use.autocorr=T #Adjust Probabilities due to Autocorrelation
diff.sd<-.025 #Width of acceleration distribution
loess.sd<-.025 #Width of loess distance distribution
probs.auto.sd<-.1 #Width of autocorrelation distribution
loess.span<-30 #Days for loess curve fit (smaller= tighter fit)
#weight of loess curve probs
#(>1 means more important than acceleration of price change)
loess.weight<-1
####Put daily data into format####
setwd(data.dir)
daily<-read.csv("goxUSDdaily.csv", header=T)
gox<-matrix(ncol=7,nrow=nrow(daily))
gox[,7]<-daily[,5]
gox2<-gox
##
###Live Plots function
plot.heatmaps.etc<-function(spot.color="Black"){
par(mfrow=c(2,5))
plot(p, type="l", lwd=3, xlim=c(length(p)-100, length(p)), ylab="USD/BTC",
ylim=c(0,max(p[(length(p)-100): length(p),])),
main=c("MtGox USD/BTC",round(tail(p,1),1)))
points(length(p),tail(p,1), col="Red", pch=16, cex=1.5)
lines(smthp, col="Blue", lwd=2)
lines(seq(length(p)-100, length(p), length=length(p)),
.8*min(p[(length(p)-100): length(p),])*p/max(p))
abline(h=max(.8*min(p[(length(p)-100): length(p),])*p/max(p)))
text(x=length(p)-80,y=.9*max(.8*min(p[(length(p)-100): length(p),])*p/max(p)),
labels=paste("Max=",round(max(p),1)), cex=.7)
plot(100*d2, type="l", xlim=c(length(p)-102, length(p)-2), ylab="% Price Change",
main=c("% Price Change",round(100*tail(d2,1),1)))
points(length(d2),100*tail(d2,1), col="Red", pch=16, cex=1.5)
abline(h=0)
par(mar=c(5.1,4.1,4.1,3.1))
img.matrix<-matrix(nrow=length(seq(min(diff(d)),max(diff(d)),by=.001)),
ncol=length(seq(min(d[-nrow(d)]),max(d[-nrow(d)]),by=.001))
)
diffd.image<-(diff(d)-min(diff(d)))/(max(diff(d))-min(diff(d)))
d.image<-(d[-nrow(d)]-min(d[-nrow(d)]))/(max(d[-nrow(d)])-min(d[-nrow(d)]))
for(i in 1:nrow(img.matrix)){
img.matrix[i,]<-dnorm(seq(min(diff(d)),max(diff(d)),by=.001),tail(diff(d2),1)[1],diff.sd)[i]
}
img.matrix<-img.matrix/max(img.matrix)
if(use.autocorr==T){
for(i in 1:ncol(img.matrix)){
img.matrix[,i]<-img.matrix[,i]+dnorm(seq(min(d[-nrow(d)]),max(d[-nrow(d)]),by=.001),delta,probs.auto.sd)[i]/
max(dnorm(seq(min(d[-nrow(d)]),max(d[-nrow(d)]),by=.001),delta,.1))
}
}
image.plot(img.matrix/max(img.matrix), axes=F, xlab="Rate of % Change", ylab="Lag-1 % Change")
points(diffd.image,d.image, pch=21, col="Black", bg="Grey")
#cur.point.x<-diffd.image[which(abs(tail(d2,1)[1]-d)==min(abs(tail(d2,1)[1]-d)))]
cur.point.x<-diffd.image[which(d==delta)]
points(cur.point.x,
(delta-min(d))/(max(d)-min(d)),
col="Grey", bg=spot.color, pch=21, cex=2)
title(main=c("Rate of Change vs", "% Change in Price"))
axis(side=1, at=seq(min(diffd.image),max(diffd.image),length=6),
100*round(seq(min(diff(d)),max(diff(d)),length=6),2)
)
axis(side=2, at=seq(min(d.image),max(d.image),length=6),
100*round(seq(min(d[-nrow(d)]),max(d[-nrow(d)]),length=6),2), las=2
)
img.matrix<-matrix(nrow=length(seq(min(p.smth.diff1),max(p.smth.diff1),by=.001)),
ncol=length(seq(min(d),max(d),by=.001))
)
p.smth.diff1.image<-(p.smth.diff1-min(p.smth.diff1))/(max(p.smth.diff1)-min(p.smth.diff1))
d.image<-(d-min(d))/(max(d)-min(d))
for(i in 1:nrow(img.matrix)){
img.matrix[i,]<-dnorm(seq(min(p.smth.diff1),max(p.smth.diff1),by=.001),
tail(p.smth.diff2,1),loess.sd)[i]
}
if(length(which(img.matrix>0))<1){
for( i in 1:length(img.matrix)){
img.matrix[i]<-1/length(img.matrix)
}
}
img.matrix<-img.matrix/max(img.matrix, na.rm=T)
if(use.autocorr==T){
for(i in 1:ncol(img.matrix)){
img.matrix[,i]<-img.matrix[,i]+dnorm(seq(min(d[-nrow(d)]),max(d[-nrow(d)]),by=.001),delta,probs.auto.sd)[i]/
max(dnorm(seq(min(d[-nrow(d)]),max(d[-nrow(d)]),by=.001),delta,.1))
}
}
image.plot(img.matrix/max(img.matrix), axes=F, xlab="% Distance to Loess Fit", ylab="Lag-1 % Change")
points(p.smth.diff1.image,d.image, pch=21, col="Black", bg="Grey")
cur.point.x<-p.smth.diff1.image[which(d==delta)]
points(cur.point.x,
(delta-min(d))/(max(d)-min(d)),
col="Grey", bg=spot.color, pch=21, cex=2)
title(main=c("% Distance to Loess Fit vs", "% Change in Price"))
axis(side=1, at=seq(min(p.smth.diff1.image),max(p.smth.diff1.image),length=6),
100*round(seq(min(p.smth.diff1),max(p.smth.diff1),length=6),2)
)
axis(side=2, at=seq(min(d.image),max(d.image),length=6),
100*round(seq(min(d[-nrow(d)]),max(d[-nrow(d)]),length=6),2), las=2
)
par(mar=c(5.1,4.1,4.1,2.1))
##############################
if(use.autocorr==T){
probs.init2<-probs.init
probs.init2[-which(d==delta)]<-probs.init[-which(d==delta)]*probs.autocorr
probs.init2[which(d==delta)]<-1/length(probs.init2)
}
normalizer<-1/(diff(range(d))/512)
y.max<-max(density(d, weights=probs.init2/sum(probs.init2))$y, normalizer*(probs.init2/sum(probs.init2)))
plot(d,normalizer*probs.init2/sum(probs.init2), ylim=c(0,y.max), xaxt="n",
xlab="Percent Price Change", ylab="Probability Density",
main="Autocorrelation Weights")
lines(density(d, weights=probs.init2/sum(probs.init2)), col="Red", lwd=3)
#points(tail(d2,1),normalizer*probs.init[which(d==delta)[1]]/sum(probs.init2),
#col="Light Blue", pch=16, cex=1.5
#)
axis(side=1, at=seq(min(d),max(d),length=6),
100*round(seq(min(d),max(d),length=6),2)
)
plot(100*diff(d2), type="l", xlim=c(length(p)-102, length(p)-2), ylab="Rate of Price Change",
main=c("Rate of % Price Change", round(tail(100*diff(d2),1),1)))
points(length(diff(d2)),100*tail(diff(d2),1), col="Red", pch=16, cex=1.5)
abline(h=0)
plot(100*p.smth.diff2, type="l", xlim=c(length(p)-101, length(p)-1), ylab="% Difference From Loess Fit",
main=c("% Difference from Loess Fit",round(tail(100*p.smth.diff2,1),1)))
abline(h=0)
points(length(p.smth.diff2),100*tail(p.smth.diff2,1), col="Red", pch=16, cex=1.5)
y.max<-max(density(d, weights=probs1/sum(probs1))$y, normalizer*(probs1/sum(probs1)))
plot(d,normalizer*probs1/sum(probs1), ylim=c(0,y.max), xaxt="n",
xlab="Percent Price Change", ylab="Probability Density",
main=c("Weighted due to", "Rate of Price Change"))
lines(density(d, weights=probs1/sum(probs1)), col="Red", lwd=3)
axis(side=1, at=seq(min(d),max(d),length=6),
100*round(seq(min(d),max(d),length=6),2)
)
#points(tail(d2,1),normalizer*probs.init[which(d==delta)[1]]/sum(probs.init),
#col="Light Blue", pch=16, cex=1.5
#)
y.max<-max(density(d, weights=probs2/sum(probs2))$y, normalizer*(probs2/sum(probs2)))
plot(d,normalizer*probs2/sum(probs2), ylim=c(0,y.max), xaxt="n",
xlab="Percent Price Change", ylab="Probability Density",
main=c("Weighted due to", "Distance from Loess Fit"))
lines(density(d, weights=probs2/sum(probs2)), col="Red", lwd=3)
axis(side=1, at=seq(min(d),max(d),length=6),
100*round(seq(min(d),max(d),length=6),2)
)
#points(tail(d2,1),normalizer*probs.init[which(d==delta)[1]]/sum(probs.init),
#col="Light Blue", pch=16, cex=1.5
#)
y.max<-max(density(d, weights=probs/sum(probs))$y, normalizer*(probs/sum(probs)))
plot(d,normalizer*probs/sum(probs), ylim=c(0,y.max), xaxt="n",
xlab="Percent Price Change", ylab="Probability Density",
main="Final Weights")
lines(density(d, weights=probs/sum(probs)), col="Red", lwd=3)
axis(side=1, at=seq(min(d),max(d),length=6),
100*round(seq(min(d),max(d),length=6),2)
)
#points(tail(d2,1),normalizer*probs.init[which(d==delta)[1]]/sum(probs.init),
#col="Light Blue", pch=16, cex=1.5
#)
#points(d,normalizer*probs4/sum(probs4), col="Grey", pch=16)
#lines(density(d, weights=probs4/sum(probs4)), col="Grey", lwd=3)
}
####End Live plots function
###Highest Density Interval Calculator Function
get.HDI<-function(sampleVec,credMass){
sortedPts = sort( sampleVec )
ciIdxInc = floor( credMass * length( sortedPts ) )
nCIs = length( sortedPts ) - ciIdxInc
ciWidth = rep( 0 , nCIs )
for ( i in 1:nCIs ) {
ciWidth[ i ] = sortedPts[ i + ciIdxInc ] - sortedPts[ i ]
}
HDImin = sortedPts[ which.min( ciWidth ) ]
HDImax = sortedPts[ which.min( ciWidth ) + ciIdxInc ]
HDIlim = c( HDImin , HDImax, credMass )
return( HDIlim )
}
###End HDI function
###Initialize variables
start.days<-seq(days.before.today.start,
0, by=-pred.days)
if(tail(start.days,1)!=0){
start.days<-c(start.days,0)
}
p.outs=NULL
means2.out=NULL
medians2.out=NULL
HDIs95.out=NULL
HDIs67.out=NULL
if(live.plots==T){
setwd(movie.dir)
#if(getwd()==movie.dir){
#do.call(file.remove, list(list.files(getwd()))
#}
png(filename = "GoxRchangePlot%03d.png", res=72,
width = 1600, height = 900, units = "px", pointsize = 24)
}
if(start.days[1]>0){back.testing=T}else{back.testing=F}
for(s in 1:length(start.days)){
start.day<-start.days[s]
gox<-gox2[1:(nrow(gox2)-start.day),]
p<-as.matrix(gox[,7]) #Get gox daily prices from imported data
d<-diff(p)/p[-nrow(p)] #Calculate day-to-day change in price
d<-as.matrix(rnorm(length(d),d,0.001))
autocorrs<-acf(d, plot=F) #Calculate autocorrelation in diffs
autocorrs<-autocorrs$acf[-1] #Remove lag-0 data from autocorrelation
probs<-rep(1,length(d))
probs.init<-probs
prior.probs<-probs.init
#prior.probs[which(d<0)]<-probs.init[which(d<0)]*2
t<-(1:length(p))
smthp<-loess(p~t,degree=1, span=loess.span/length(p),evaluation=length(p),
family="gaussian", control = loess.control(surface = "direct"))$fitted
p.smth.diff1<-(p-smthp)/smthp
p.smth.diff1<-p.smth.diff1[-length(p.smth.diff1)]
####################
#Allocate Memory for output matrix and initialize progress bar
p.out=matrix(nrow=nrow(p)+pred.days, ncol=iters)
pb<-txtProgressBar(min = 0, max = iters,
initial = 0, char = "=",width = 60, style = 3,)
#
#Run simulation
for(j in 1:iters){
p<-as.matrix(gox[,7])
t<-(1:length(p))
smthp<-loess(p~t,degree=1, span=loess.span/length(p),evaluation=length(p),
family="gaussian", control = loess.control(surface = "direct"))$fitted
p.smth.diff1<-(p-smthp)/smthp
p.smth.diff1<-p.smth.diff1[-length(p.smth.diff1)]
delta<-tail(d,1)[1]
for(i in 1:pred.days){
probs1<-probs.init
probs2<-probs.init
probs.autocorr<-probs.init
probs<-probs.init
d2<-diff(p)/p[-nrow(p)]
dtemp<-cbind(d[-nrow(d)],diff(d))
dtemp2<-cbind(order(dtemp[,2]),dtemp[order(dtemp[,2]),])
dtemp2<-cbind(dtemp2,dnorm(dtemp2[,3],tail(diff(d2),1),diff.sd))
probs1[-nrow(d)]<-dtemp2[order(dtemp2[,1]),4]
probs1<-probs1/sum(probs1)
probs1[nrow(probs1)]<-1/length(probs)
t<-(1:length(p))
smthp<-loess(p~t,degree=1, span=loess.span/length(p),evaluation=length(p),
family="gaussian", control = loess.control(surface = "direct"))$fitted
p.smth.diff2<-(p-smthp)/smthp
dtemp<-cbind(d,p.smth.diff1)
dtemp3<-cbind(order(dtemp[,2]),dtemp[order(dtemp[,2]),])
dtemp3<-cbind(dtemp3,dnorm(dtemp3[,3],tail(p.smth.diff2,1),loess.sd))
probs2<-dtemp3[order(dtemp3[,1]),4]
probs2<-probs2/sum(probs2)
probs<-probs1+loess.weight*probs2
if(use.autocorr==T){
probs.autocorr<-dnorm(d[-which(d==delta)],delta,probs.auto.sd)
probs.autocorr<-probs.autocorr/sum(probs.autocorr)
probs[-which(d==delta)]<-probs[-which(d==delta)]+ probs.autocorr
}
probs<-probs*prior.probs
if(live.plots==T){
if(j==1 && s==1){
plot.heatmaps.etc(spot.color="Black")
}
}
if(length(which(probs>0))<1){probs<-probs.init}
delta<-sample(d,1, prob=probs)
if(live.plots==T){
if(j==1 && s==1){
plot.heatmaps.etc(spot.color="Blue")
}
}
new<-as.numeric(tail(p,1)+delta*tail(p,1))
if(new<0){new=0}
p<-rbind(p,new[1])
}
setTxtProgressBar(pb, j)
p.out[,j]<-p
if(live.plots==T && s==1 && j==1){dev.off()
setwd(data.dir)
}
}
close(pb)
print(paste(s, "of", length(start.days), "Complete"))
#End Simulation
#Calculate means, medians, and hdis at each timepoint
means=matrix(nrow=nrow(p.out),ncol=1)
medians=matrix(nrow=nrow(p.out),ncol=1)
HDIs95=matrix(nrow=nrow(p.out),ncol=2)
HDIs67=matrix(nrow=nrow(p.out),ncol=2)
for(i in 1:nrow(p.out)){
means[i,]<-mean(p.out[i,])
medians[i,]<-median(p.out[i,])
HDIs95[i,]<-cbind(get.HDI(p.out[i,], .95)[1],get.HDI(p.out[i,], .95)[2])
HDIs67[i,]<-cbind(get.HDI(p.out[i,], .67)[1],get.HDI(p.out[i,], .67)[2])
}
#
p.out2<-p.out[-(1:nrow(gox)),]
means2<-means[-(1:nrow(gox)),]
medians2<-medians[-(1:nrow(gox)),]
HDIs95<-HDIs95[-(1:nrow(gox)),]
HDIs67<-HDIs67[-(1:nrow(gox)),]
p.outs<-rbind(p.outs,p.out2)
means2.out<-c(means2.out,means2)
medians2.out<-c(medians2.out,medians2)
HDIs95.out<-rbind(HDIs95.out,HDIs95)
HDIs67.out<-rbind(HDIs67.out,HDIs67)
}
gox<-gox2[1:(nrow(gox2)-start.days[1]),]
means<-as.matrix(c(gox[,7],means2.out))
medians<-as.matrix(c(gox[,7],medians2.out))
HDIs95<-rbind(cbind(gox[,7],gox[,7]),HDIs95.out)
HDIs67<-rbind(cbind(gox[,7],gox[,7]),HDIs67.out)
goxtemp=matrix(nrow=nrow(gox),ncol=ncol(p.out))
goxtemp[,1:ncol(goxtemp)]<-gox[,7]
p.out<-rbind(goxtemp,p.outs)
#############################
####Plot predictions and posterior distributions
if(make.pngs==T){
png(filename = "GoxPredict%03d.png", res=150,
width = 1920, height = 1080, units = "px", pointsize = 12)
}else{
dev.new()
}
layout(matrix(c(1,1,2,3,4,5,6,7,8),nrow=3, ncol=3, byrow=T))
if(
diff(range(c(gox[max(1,nrow(gox)-3*pred.days):nrow(gox),7],
means[nrow(gox):nrow(means),])))>300
){
plot(0,0.1,
#ylim=c(0.1,max(HDIs95[,2])), xaxt="n",
ylim=c(0.1,max(HDIs95[,2])), xaxt="n", log="y",
xlim=c(max(0,nrow(gox)-3*pred.days), nrow(p.out)), type="n",
#xlab=paste("Days From 'Today' (", start.days[1], " Days Ago)", sep=""),
xlab=paste("Days From Today"),
ylab="USD/BTC", main=c("MtGox USD/BTC Predictions",
paste("Iterations:", iters)
), sub=paste("Use AutoCorrelation Data=", use.autocorr))
axis(side=1,
at=floor(seq(max(0,nrow(gox)-3*pred.days), nrow(p.out), length=6)),
labels=floor(seq(max(0,nrow(gox)-3*pred.days), nrow(p.out), length=6)-nrow(gox2)
))
}else{
plot(0,0.1,
ylim=c(0.01,max(HDIs95[,2])), xaxt="n",
#ylim=c(0.01,max(HDIs95[,2])), xaxt="n", log="y",
xlim=c(max(0,nrow(gox)-3*pred.days), nrow(p.out)), type="n",
#xlab=paste("Days From 'Today' (", start.days[1], " Days Ago)", sep=""),
xlab=paste("Days From Today"),
ylab="USD/BTC", main=c("MtGox USD/BTC Predictions",
paste("Iterations:", iters)
), sub=paste("Use AutoCorrelation Data=", use.autocorr))
axis(side=1,
at=floor(seq(max(0,nrow(gox)-3*pred.days), nrow(p.out), length=6)),
labels=floor(seq(max(0,nrow(gox)-3*pred.days), nrow(p.out), length=6)-nrow(gox2)
))
}
for(i in seq(1,ncol(p.out),length=1000)){
lines(p.out[,i],type="l", col=rainbow(iters, alpha=1/min(iters,10))[i])
}
legend("topleft", c("Mean","Median","95% HDI","67% HDI"),
col=c("Black","Dark Grey","Black","Light Grey"), lty=c(1,1,2,2),
bty="n", lwd=2, ncol=2, text.width=20, cex=.7
)
abline(v=nrow(gox))
abline(v=nrow(gox2))
points(nrow(gox)+floor(seq(1, start.days[1]+pred.days, length=6)),rep(.1,6), pch=2)
lines(nrow(as.matrix(gox[,7])): nrow(p.out),
means[nrow(as.matrix(gox[,7])): nrow(p.out),], col="Black", lwd=3)
lines(nrow(as.matrix(gox[,7])): nrow(p.out),
medians[nrow(as.matrix(gox[,7])): nrow(p.out),], col="Light Grey", lwd=3)
lines(nrow(as.matrix(gox[,7])): nrow(p.out),
HDIs95[nrow(as.matrix(gox[,7])): nrow(p.out),1], col="Black", lwd=3, lty=2)
lines(nrow(as.matrix(gox[,7])): nrow(p.out),
HDIs95[nrow(as.matrix(gox[,7])): nrow(p.out),2], col="Black", lwd=3, lty=2)
lines(nrow(as.matrix(gox[,7])): nrow(p.out),
HDIs67[nrow(as.matrix(gox[,7])): nrow(p.out),1], col="Light Grey", lwd=3, lty=2)
lines(nrow(as.matrix(gox[,7])): nrow(p.out),
HDIs67[nrow(as.matrix(gox[,7])): nrow(p.out),2], col="Light Grey", lwd=3, lty=2)
lines(as.matrix(gox2[,7]), lwd=3)
lines(as.matrix(gox2[,7]), lwd=1, col="Blue")
d.temp<-d
hist(100*d.temp, freq=F, col="Green", ylim=c(0,max(density(100*d)$y)),
breaks=seq(min(100*d.temp),1.1*max(100*d.temp),
by=.025*diff(range(100*d.temp))),
xlab="Percent Change",
main=c("Day-to-day % Change",
paste("95% HDI:",round(get.HDI(100*d.temp,.95)[1],1),
"-",round(get.HDI(100*d.temp,.95)[2],1)),
paste("67% HDI:",round(get.HDI(100*d.temp,.67)[1],1),
"-",round(get.HDI(100*d.temp,.67)[2],1))
),
sub=paste("Mean=",round(mean(100*d.temp),2), "Median=", round(median(100*d.temp),2))
)
lines(density(100*d.temp), lwd=3)
segments(get.HDI(100*d.temp,.95)[1],0,get.HDI(100*d.temp,.95)[2],0, lwd=5)
segments(get.HDI(100*d.temp,.67)[1],0,get.HDI(100*d.temp,.67)[2],0, col="Light Grey", lwd=3)
for(i in floor(seq(1, start.days[1]+pred.days, length=6))){
y.max<-max(density(p.out[nrow(as.matrix(gox[,7]))+i,])$y)
x.lim<-c(min(p.out[nrow(as.matrix(gox[,7]))+i,]),
HDIs95[nrow(as.matrix(gox[,7]))+i,][2])
hist(p.out[nrow(as.matrix(gox[,7]))+i,], freq=F,
xlab="USD/BTC", col="Grey",ylim=c(0,y.max),
xlim=x.lim,
breaks=seq(min(p.out[nrow(as.matrix(gox[,7]))+i,]),
1.1*max(p.out[nrow(as.matrix(gox[,7]))+i,]),
by=.05*diff(x.lim)),
sub=paste("Mean=", round(means[nrow(as.matrix(gox[,7]))+i,],2),
"Median=", round(medians[nrow(as.matrix(gox[,7]))+i,],2)
), main=""
)
title(main=c(paste(-(start.days[1]-i), "Days From Today"),
paste("Range:",round(min(p.out[nrow(as.matrix(gox[,7]))+i,]),1),"-",
round(max(p.out[nrow(as.matrix(gox[,7]))+i,]),1)),
paste("95% HDI:", round(HDIs95[nrow(as.matrix(gox[,7]))+i,][1],1), "-",
round(HDIs95[nrow(as.matrix(gox[,7]))+i,][2],1)),
paste("67% HDI:", round(HDIs67[nrow(as.matrix(gox[,7]))+i,][1],1), "-",
round(HDIs67[nrow(as.matrix(gox[,7]))+i,][2],1))
),
cex.main=1)
lines(density(p.out[nrow(as.matrix(gox[,7]))+i,], n=1024),
col="Black", lwd=3
)
if(back.testing==T){
if((nrow(as.matrix(gox[,7]))+i) if(is.finite(gox2[nrow(as.matrix(gox[,7]))+i,7])){
if(HDIs95[nrow(as.matrix(gox[,7]))+i,][1] HDIs95[nrow(as.matrix(gox[,7]))+i,][2]>gox2[nrow(as.matrix(gox[,7]))+i,7]){
y.dens<-density(p.out[nrow(as.matrix(gox[,7]))+i,], n=1024)$y[which(density(p.out[nrow(as.matrix(gox[,7]))+i,], n=1024)$x>
gox2[nrow(as.matrix(gox[,7]))+i,7])[1]]
points(gox2[nrow(as.matrix(gox[,7]))+i,7],y.dens, pch=16, col="Red", cex=1.5)
text(gox2[nrow(as.matrix(gox[,7]))+i,7],y.dens,gox2[nrow(as.matrix(gox[,7]))+i,7], pos=4,
offset=1, col="Blue", font=2)
}else{
if(gox2[nrow(as.matrix(gox[,7]))+i,7] text(x.lim[1],y.max,gox2[nrow(as.matrix(gox[,7]))+i,7], pos=4,
offset=1, col="Red", font=2)
}else{
text(.9*x.lim[2],y.max,gox2[nrow(as.matrix(gox[,7]))+i,7], pos=4,
offset=1, col="Red", font=2)
}
}
}
}
}
segments(HDIs95[nrow(as.matrix(gox[,7]))+i,][1], 0,
HDIs95[nrow(as.matrix(gox[,7]))+i,][2],0, lwd=5)
segments(HDIs67[nrow(as.matrix(gox[,7]))+i,][1], 0,
HDIs67[nrow(as.matrix(gox[,7]))+i,][2],0, col="Light Grey", lwd=3)
}
if(make.pngs==T){dev.off()}
#
#Plot example traces
if(make.pngs==T){
png(filename = "GoxExamples%03d.png", res=150,
width = 1920, height = 1080, units = "px", pointsize = 12)
}else{
dev.new()
}
par(mfrow=c(1,1))
examps<-floor(runif(num.examps,1,ncol(p.out)))
plot(0,0.1,
ylim=c(0.001,max(p.out[,examps],as.matrix(gox2[,7]))), xaxt="n",
xlim=c(max(0,nrow(gox)-3*pred.days), nrow(p.out)), type="n",
xlab="Days From Today",ylab="USD/BTC",
main=c("MtGox USD/BTC Predictions", paste(num.examps,"Example Scenarios"))
)
for(i in examps){
lines(p.out[,i],type="l", col=rainbow(num.examps, alpha=1)[which(examps==i)], lwd=2)
}
axis(side=1,
at=floor(seq(max(0,nrow(gox)-3*pred.days), nrow(p.out), length=6)),
labels=floor(seq(max(0,nrow(gox)-3*pred.days), nrow(p.out), length=6)-nrow(gox)
))
if(back.testing==T){
lines(as.matrix(gox2[,7]), lwd=5)
}else{
lines(as.matrix(gox[,7]), lwd=5)
}
abline(v=nrow(gox))
if(make.pngs==T){dev.off()}
#
#graphics.off()
Update the prediction one step into the future and repeat 7
