We shall denote the integral with the letter $I$, i.e. :

$$I:=\int_0^\infty \frac{\mathrm{d}x}{x} \int_0^x \frac{\mathrm{d}y}{y} \int_0^y \frac{\mathrm{d}z}{z} [\sin{x}+\sin{(x-y)}-\sin{(x-z)}-\sin{(x-y+z)}]$$

The key point is to insert a damping exponential $e^{-\alpha z}$ (where $\alpha \geq 0$) into the integrand, this makes the integral directly dependent on the $\alpha$ factor, then apply the Leibnitz formula for differentatiaon under the integral sign, finally integrate again recalling $I(\infty)=0$ and $I(0)=I$. Albeit one can view on this procedure, respectively it can be done, using replacement $$\frac{1}{z}=\int_0^\infty \! e^{-\alpha z} \; {\mathrm{d}\alpha}$$

Then, the change of the order of integration would imply this reformulation to appear :

$$I=\int_0^\infty \! {\mathrm{d}\alpha} \int_0^\infty \frac{\mathrm{d}x}{x} \int_0^x \frac{\mathrm{d}y}{y} \int_0^y {\mathrm{d}z} \, [\sin{x}+\sin{(x-y)}-\sin{(x-z)}-\sin{(x-y+z)}] \, e^{-\alpha z}$$

The inner integral is now directrly evaluable, holding the line :

$$I=\int_0^\infty \! {\mathrm{d}\alpha} \! \int_0^\infty \! \frac{\mathrm{d}x}{x} \int_0^x \frac{\mathrm{d}y}{y} \! \frac{1-e^{-\alpha y}}{\alpha(1+\alpha^2)} (\sin{x}+\sin{(x-y)}) +\frac{1+e^{-\alpha y}}{1+\alpha^2}(\cos{x}-\cos{(x-y)}) $$

As you can fine in this article (http://arxiv.org/pdf/1201.1975v1.pdf), integral similar to one contained in $I$ is computed there, namely :

$$\zeta{(2)} = \int_0^\infty \! \frac{\mathrm{d}x}{x} \int_0^x \frac{\mathrm{d}y}{y} (\cos{(x-y)}-\cos{x}) $$

Althought, from representation of $I$ we are able to derive this result (denoted $\tilde{I}$ in the article) regardless of using the final one as a fact or from a different source. The actual approach is however very similar as the one described

in the article : Just consider a parametric integral $\tilde{I}(\alpha)$ defined as:

$$\tilde{I}(\alpha) := \int_0^\infty \! \frac{\mathrm{d}x}{x} \int_0^x \frac{1+e^{-\alpha y}}{y} (\cos{(x-y)}-\cos{x}) \, \mathrm{d}y $$

Evidently : $\tilde{I}(0)=2\tilde{I} \wedge \tilde{I}(\infty)=\tilde{I}$, using Fundamental Theorem of Calculus we can write $$\tilde{I}=-\left( \tilde{I}-2\tilde{I} \right)=-\left( \tilde{I}(\infty)-\tilde{I}(0) \right)=-\int_0^\infty \tilde{I'}(\alpha) \, \mathrm{d}\alpha$$

Then differentiation under the integral sign and some tricky integration will lead us to the value $\zeta{(2)}$ as well.

Nevertheless, we will now continue in evaluation of $I$, the effort done on the discussion of the $\zeta{(2)}$ integral is not meaningless at all, because now we will separate this result from our integral representation of $I$, because of the relation : $$1+e^{-\alpha y}=2-(1-e^{-\alpha y})$$

we now split the $I$-integral into two pieces :

$$I=J-2\int_0^\infty \mathrm{d}\alpha \frac{\zeta{(2)}}{1+\alpha^2} =J-\frac{\pi^3}{6}$$

where

$$J:=\int_0^\infty \! {\mathrm{d}\alpha} \! \int_0^\infty \! \frac{\mathrm{d}x}{x} \int_0^x \frac{\mathrm{d}y}{y} \! \frac{1-e^{-\alpha y}}{\alpha(1+\alpha^2)} (\sin{x}+\sin{(x-y)})-\frac{1-e^{-\alpha y}}{1+\alpha^2}(\cos{x}-\cos{(x-y)})$$

We are doing this for preparation to make an unified replacement (another), namely :

$$\frac{1-e^{-\alpha y}}{y}=\int_0^\alpha \! e^{-\beta y} \; {\mathrm{d}\beta}$$

After change of the order of integration in $\beta$ :

$$J=\int_0^\infty \! \frac{\mathrm{d}\alpha}{\alpha(1+\alpha^2)} \! \int_0^\alpha \! {\mathrm{d}\beta} \! \int_0^\infty \! \frac{\mathrm{d}x}{x} {\int_0^x} \mathrm{d}y \,[\sin{x}+\sin{(x-y)}-\alpha (\cos{x}-\cos{(x-y)})] e^{-\beta y}$$

Again, by direct integration of inner integral :

$$J\!=\!\int_0^\infty \!\!\!\! \frac{\mathrm{d}\alpha}{\alpha(1+\alpha^2)} \!\! \int_0^\alpha \!\! {\mathrm{d}\beta} \! \int_0^\infty \!\! \mathrm{d}x \frac{\sin{x}\!-\!\alpha\cos{x}}{\beta}\frac{1\!-\!e^{-\beta y}}{x}\!+\!\frac{\alpha \beta-\!1}{1\!+\!\beta^2}\frac{\cos{x}\!-\!e^{-\beta y}}{x}\!+\!\frac{\alpha\!+\!\beta}{1\!+\!\beta^2}\frac{\sin{x}}{x}$$

Resulting integrals are somewhat elementary (i.e. Dirichlet integral etc.), however we can use that dumb replacement rule again using

$$\frac{1}{x}=\int_0^\infty \! e^{-\gamma x} \; {\mathrm{d}\gamma}$$

The other replacement however could also be done, i.e.

$$\frac{1-e^{-\beta x}}{x}=\int_0^\beta \! e^{-\gamma x} \; {\mathrm{d}\gamma}$$

Define $$ \begin{split} S(\beta) & :=\int_0^\infty \mathrm{d}x \frac{1-e^{-\beta y}}{x}\sin{x} \\ C(\beta) & :=\int_0^\infty \mathrm{d}x \frac{1-e^{-\beta y}}{x}\cos{x} \\ L(\beta) & :=\int_0^\infty \mathrm{d}x \frac{\cos{x}-e^{-\beta y}}{x} \\ D & :=\int_0^\infty \mathrm{d}x \frac{\sin{x}}{x} \end{split}$$

By definition then

$$ J = \int_0^\infty \frac{\mathrm{d}\alpha}{\alpha(1+\alpha^2)} \int_0^\alpha {\mathrm{d}\beta} \; \frac{1}{\beta}S(\beta)-\frac{\alpha}{\beta} C(\beta)+\frac{\alpha \beta-1}{1+\beta^2} L(\beta)+\!\frac{\alpha+\beta}{1+\beta^2} D $$

Now on the evaluation - for the first two integrals $S$ and $C$ the quickest way to compute them is via second type of replacement, The third an fourth one are done (due to the convergence issues) by the "dumb" replacement. For each integral we have: $$S(\beta)=\int_0^\beta {\mathrm{d}\gamma} \int_0^\infty \mathrm{d}x \sin{x} e^{-\gamma x} = \int_0^\beta {\mathrm{d}\gamma} \frac{1}{1+\gamma^2}=\arctan{\beta} \\ \\ C(\beta)=\int_0^\beta {\mathrm{d}\gamma} \int_0^\infty \mathrm{d}x \cos{x} e^{-\gamma x} = \int_0^\beta {\mathrm{d}\gamma} \frac{\gamma}{1+\gamma^2}=\frac{1}{2} \ln{(1+\beta^2)} \\ L(\beta)=\int_0^\infty {\mathrm{d}\gamma} \int_0^\infty \mathrm{d}x (\cos{x}\!-\!e^{-\beta x}) e^{-\gamma x} = \int_0^\infty {\mathrm{d}\gamma} \frac{\gamma}{1\!+\!\gamma^2}\!-\!\frac{1}{\beta\!+\!\gamma} = \ln{\frac{\sqrt{1 \!+\!\gamma^2}}{\beta \! + \! \gamma}} {\Big|}_0^\infty = \ln{\beta} \\ D=\int_0^\infty {\mathrm{d}\gamma} \int_0^\infty \mathrm{d}x \sin{x} e^{-\gamma x} = \int_0^\infty \frac{\mathrm{d}\gamma}{1+\gamma^2} = \frac{\pi}{2} $$

Substituting these results to the $J$ integral :

$$ J = \int_0^\infty \frac{\mathrm{d}\alpha}{\alpha(1+\alpha^2)} \int_0^\alpha {\mathrm{d}\beta} \; \frac{\arctan{\beta}}{\beta}-\frac{\alpha}{2\beta} \ln{(1+\beta^2)}+\frac{\alpha \beta-1}{1+\beta^2} \ln{\beta}+\!\frac{\alpha+\beta}{1+\beta^2} \frac{\pi}{2} $$

Next step will turn out to be really important. Just consider how far we are, we ended up with a double integral with its domain of integration - infinite triangle, i.e. half of the first quadrant. Let us perform a variable change in such manner it transforms to rectangle. There are many candidates (e.g. polar one), however we chose this particular :

$$\alpha=t \\ \beta=tk$$

With Jacobian of transformation equal to $J=t$, domain itself now consists of a semiinfinite strip $$(t,k)\in\{(0,\infty)\times(0,1)\}$$

After this transformation :

$$ J \!=\! \int_0^1 \!\mathrm{d}k \int_0^\infty \!\mathrm{d}t \; \frac{\arctan{kt}}{kt(1+t^2)}-\frac{\ln{(1\!+\!k^2t^2)}}{2k(1\!+\!t^2)}\!-\!\frac{1-kt^2}{(1\!+\!t^2)(1\!+\!k^2t^2)} \ln{(kt)}\!+\!\!\frac{t(1+k)}{(1\!+\!t^2)(1\!+\!k^2t^2)} \frac{\pi}{2} $$

Defining ($0<s<4$) :

$$ \begin{split} J_1(k) & :=\int_0^\infty \frac{\arctan{kt}}{t(1+t^2)} \mathrm{d}t \\ J_2(k) & :=\int_0^\infty \frac{\ln{(1\!+\!k^2t^2)}}{2(1+t^2)} \mathrm{d}t \\ P_s(k) & :=\int_0^\infty \frac{t^{s-1}\,\mathrm{d}t}{(1\!+\!t^2)(1\!+\!k^2t^2)} \\ L_s(k) & :=\int_0^\infty \frac{t^{s-1}\ln{t} \, \mathrm{d}t}{(1+t^2)(1+k^2t^2)} \end{split} $$

Then by direct definition :

$$ J \!=\! \int_0^1 \!\mathrm{d}k \; \frac{J_1(k)}{k}-\frac{J_2(k)}{k}\!-P_1(k) \ln{k}+ P_3(k) k \ln{k} - L_1(k) + L_3(k) k + \frac{\pi}{2} (1+k)P_2(k) $$

Let us compute the $P_s(k)$ value only (for arbitrary $s$), note that :

$$\frac{1}{(1+t^2)(1+k^2t^2)}=\frac{1}{1-k^2}\frac{1}{1+t^2}-\frac{k^2}{1-k^2}\frac{1}{1+k^2t^2}$$

Substituting that decomposition to the definition of $P_s(k)$ and by inserting the substitution $tk \rightarrow t$ into the second integral we can get the $s$-dependency out of both integrals, recollecting :

$$P_s(k) =\frac{1-k^{2-s}}{1-k^2} \int_0^\infty \frac{t^{s-1}}{1+t^2} \mathrm{d}t$$

Unfortunately this representation restricts the value of $s$ greatly, now $0<s<2$

Recall the definition of the beta function, using a simple tranformation $1+t^2 \rightarrow 1/t$ and applying the reflection formula one gets:

$$ P_s(k) =\frac{1-k^{2-s}}{1-k^2} \frac{\pi}{2 \sin{(\frac{\pi s}{2})}} $$

For $s=1$ we have :

$$P_1(k) = \frac{\pi}{2} \frac{1}{1+k}$$

Using the L'Hopital's rule also extraction of $P_2(k)$ is possible :

$$P_2(k) = - \frac{\ln{k}}{1-k^2}$$

From the definition is clear that $P_1(k)+k^2P_3(k)=\int_0^\infty \! \frac{\mathrm{d}t}{1+t^2} = \pi/2$, therefore :

$$P_3(k) = \frac{\pi}{2k} \frac{1}{1+k}$$

Amazingly this result could be obtained also by evaluating the $P_s$ fromula at $s=3$

Evidently $L_s(k)=\frac{\mathrm{d}}{\mathrm{d}s}P_s(k)$, performing that differentiation on the precomputed $P_s(k)$ one gets :

$$ L_s(k) =\frac{k^{2-s} \ln{k}}{1-k^2} \frac{\pi}{2 \sin{(\frac{\pi s}{2})}} - \frac{\pi^2}{4} \frac{1-k^{2-s}}{1-k^2} \frac{\cos{(\frac{\pi s}{2})}}{\sin^2{(\frac{\pi s}{2})}} $$

Evaluation at $s=1$ gives:

$$L_1(k) = \frac{\pi}{2} \frac{k \ln{k}}{1-k^2}$$

Using the same trick as in $P_3$ one can also evalute $L_3$ (However direct substitution also could have been done)

$$L_3(k) = - \frac{\pi}{2k} \frac{\ln{k}}{1-k^2}$$

Let us now make a closer look at $J_1$ and $J_2$ integrals, as we can see, both at $k=0$ vanish, i.e. $J_1(0)=J_2(0)=0$

It seems natural to apply differentiation under the integral sign, performing that one can obtain :

$$J'_1(k) = \int_0^\infty \frac{\mathrm{d}t}{(1+t^2)(1+k^2t^2)} = P_1(k) = \frac{\pi}{2} \frac{1}{1+k} \\ J'_2(k) = k\int_0^\infty \frac{t^2 \mathrm{d}t}{(1+t^2)(1+k^2t^2)} = k P_3(k) = \frac{\pi}{2} \frac{1}{1+k} $$

Using Fundamental theorem of calculus $f(x)=f(0)+\int_0^x f'(y) \mathrm{d}y$, resulting in :

$$J_1(k)=J_2(k)= \frac{\pi}{2}\ln{(1+k)}$$

Putting all these informations ($J_s,P_s,L_s$) togehter, in $J$ we have now, due to "huge cancelerin" :

$$ J = -\pi \int_0^1 \! \frac{\ln{k}}{1-k} \, \mathrm{d}k = -\pi \mathrm{Li}_2(1-k) {\Big|}_0^1 = \frac{\pi^3}{6} $$

Tracing back to the original definition of $J$ gives the desired result :) ...