diff --git a/demos/octave/SchreiberTransferEntropyExamples/runHeartBreathRateKraskov.m b/demos/octave/SchreiberTransferEntropyExamples/runHeartBreathRateKraskov.m
new file mode 100755
index 0000000..0b1b779
--- /dev/null
+++ b/demos/octave/SchreiberTransferEntropyExamples/runHeartBreathRateKraskov.m
@@ -0,0 +1,73 @@
+% function [teHeartToBreath, teBreathToHeart] = runHeartBreathRateKraskov(knns)
+%
+% runHeartBreathRateKraskov
+% Version 1.0
+% Joseph Lizier
+% 22/1/2014
+%
+% Used to explore information transfer in the heart rate / breath rate example of Schreiber --
+%  but estimates TE using Kraskov-Grassberger estimation.
+%
+% Runtime is ~3 minutes per knn value (as there is no fast nearest neighbour search implemented yet in the toolkit)
+%
+% Inputs
+% - knns - a scalar specifying a single, or vector specifying multiple, value of K nearest neighbours to evaluate TE (Kraskov) with
+% Outputs
+% - teHeartToBreath - TE (heart -> breath) for each value of k nearest neighbours
+% - teBreathToHeart - TE (breath -> heart) for each value of k nearest neighbours
+
+
+function [teHeartToBreath, teBreathToHeart] = runHeartBreathRate(knns)
+
+	starttime = time;
+	
+	% Add utilities to the path
+	addpath('..');
+
+	% Assumes the jar is two levels up - change this if this is not the case
+	% Octave is happy to have the path added multiple times; I'm unsure if this is true for matlab
+	javaaddpath('../../../infodynamics.jar');
+
+	fprintf('TE for heart rate <-> breath rate for Kraskov estimation:\n');
+
+	data = load('../../data/SFI-heartRate_breathVol_bloodOx.txt');
+	
+	% Restrict to the samples that Schreiber mentions:
+	heartRate = data(2350:3550,:);
+	
+	% Separate the data from each column:
+	heart = data(:,1);
+	chestVol = data(:,2);
+	bloodOx = data(:,3);
+	timeSteps = length(heart);
+	
+	% Using a single conditional mutual information calculator is the least biased way to run this:
+	teCalc=javaObject('infodynamics.measures.continuous.kraskov.ConditionalMutualInfoCalculatorMultiVariateKraskov1');
+
+	for knnIndex = 1:length(knns)
+		knn = knns(knnIndex);
+		% Compute a TE value for knn nearest neighbours
+				
+		% Perform calculation for heart -> breath (lag 1)
+		teCalc.initialise(1,1,1); % Use history length 1 (Schreiber k)
+		teCalc.setProperty('k', sprintf('%d',knn));
+		teCalc.setObservations(octaveToJavaDoubleMatrix(heart(1:timeSteps-1)), ...
+					octaveToJavaDoubleMatrix(chestVol(2:timeSteps)), ...
+					octaveToJavaDoubleMatrix(chestVol(1:timeSteps-1)));
+		teHeartToBreath(knnIndex) = teCalc.computeAverageLocalOfObservations();
+		
+		% Perform calculation for breath -> heart (lag 1)
+		teCalc.initialise(1,1,1); % Use history length 1 (Schreiber k)
+		teCalc.setProperty('k', sprintf('%d',knn));
+		teCalc.setObservations(octaveToJavaDoubleMatrix(chestVol(1:timeSteps-1)), ...
+					octaveToJavaDoubleMatrix(heart(2:timeSteps)), ...
+					octaveToJavaDoubleMatrix(heart(1:timeSteps-1)));
+		teBreathToHeart(knnIndex) = teCalc.computeAverageLocalOfObservations();
+		
+		fprintf('TE(k=%d): heart->breath = %.3f, breath->heart = %.3f\n', knn, teHeartToBreath(knnIndex), teBreathToHeart(knnIndex));
+	end
+		
+	endtime = time;
+	printf("Total runtime was %.1f sec\n", endtime - starttime);
+end
+