Function RndExp(Avg As Single)
r = Rnd()
Do While r = 0# Or r = 1#
  r = Rnd()
Loop
RndExp = -Avg * Log(r)
End Function

Function GenData(Optional printIt As Integer = 0) As Integer()
Const nSyst As Integer = 4 
Const nStates = 5
Const nPeriods As Integer = 50
Dim s As String
Dim lambda As Single
Dim v As Single
Dim tau As Single
Dim nxtJump(1 To nSyst) As Single ' For each system, we keep track of the NEXT transition time, which either is before or after the next inspection
Dim t As Single ' Running time for the simulation
Dim needUpdate As Boolean
tau = Range("maxtime") ' Read from spreadsheet
Dim i As Integer, j As Integer ' i = system number, j = period number
Dim dta() As Integer
ReDim dta(1 To nSyst, 1 To nPeriods + 1)
lambda = Range("lambda")
v = Range("v")
' Initialize, and pick first transition time for each system
s = Format(1)
For i = 1 To nSyst
  dta(i, 1) = IIf(i < nStates, i, 1) ' Not all systems starts in the same state
  nxtJump(i) = RndExp(1# / (lambda * (1# + v) ^ dta(i, 1))) ' Pick the first transition time for each system
  dta(i, 2) = dta(i, 1)
  s = s & "  " & Format(dta(i, 1))
Next i
If printIt Then
Debug.Print "P# S1 S2 S3 S4"
Debug.Print "-- -- -- -- --"
Debug.Print s
End If
t = 0
For j = 2 To nPeriods
  t = t + tau
  needUpdate = True
  Do While needUpdate
    needUpdate = False
    For i = 1 To nSyst
      If nxtJump(i) < t Then
        ' A transition has taken place since last inspection
        needUpdate = True
        dta(i, j) = dta(i, j) + 1
        nxtJump(i) = nxtJump(i) + RndExp(1# / (lambda * (1# + v) ^ dta(i, 1))) ' Pick next transition time for this system
      End If
    Next i
    ' In principle there could be more than one transition for the same system since last inspection
  Loop
  s = Format(j)
  For i = 1 To nSyst
    ' Chekc for the need to repair
    If dta(i, j) > nStates Then dta(i, j) = nStates ' It does not make sense to go beyond nStates (i.e., r)
    s = s & "  " & Format(dta(i, j))
    If dta(i, j) >= nStates - 1 Then
      dta(i, j + 1) = 1
      nxtJump(i) = t + RndExp(1# / lambda) ' Pick next transition time
    Else
      dta(i, j + 1) = dta(i, j)
    End If
  Next i
  If printIt Then Debug.Print s
Next j
GenData = dta
End Function


 
P# S1 S2 S3 S4
-- -- -- -- --
1  1  2  3  4
2  3  4  4  5x
3  3  5  2  2
4  3  2  3  5
5  3  5  3  1
6  3  2  4  2
7  4  3  1  4
8  2  3  2  2
9  2  4  4  2
10  2  3  1  5
11  2  4  2  1
12  5  2  4  1
13  1  4  3  4
14  1  2  5  1
15  3  3  1  1
16  4  5  1  1
17  1  1  1  2
18  2  1  2  5
19  4  1  5  3
20  2  2  2  3
21  4  2  3  4
22  4  3  5  3
23  1  5  1  5
24  3  2  1  2
25  5  4  3  4
26  2  1  3  1
27  4  2  3  1
28  2  4  4  2
29  3  1  3  2
30  5  1  4  4
31  1  1  2  1
32  1  2  3  2
33  1  2  4  5
34  2  2  3  2
35  2  3  5  5
36  3  5  2  3
37  4  1  3  4
38  2  3  4  1
39  2  4  2  1
40  5  1  3  2
41  1  4  5  2
42  2  3  2  3
43  4  5  4  4
44  1  2  1  2
45  2  3  2  4
46  2  5  3  1
47  3  1  4  2
48  3  1  1  3
49  3  1  2  5
50  4  2  5  1