Julia Community 🟣: Matthew leung

Playing with Broadcast function

Matthew leung — Sun, 07 Aug 2022 07:10:31 +0000

One of the features in Julia is the broadcast function. It is very easy to apply the function to each element in the array with the broadcast operator ".".

However, if the function needs to take some other arguments which is not for iteration, how to do that?

Let's take an example. How to use broadcast function to calculate the moving average?

The for-loop version of the moving average is:

function ma_loop(vec,win_size)
    n = size(vec)[1]
    ma = Array{Float64}(undef, n)
    for i in 1:n
        if i >= win_size
            ma[i] = sum(vec[i-win_size+1:i])/win_size
        end
    end
    return ma[win_size:end]
end

To use broadcast function to replace the for-loop, I try to iterate the index array (1:n), and passing the Reference of the whole vector (Ref(vec)) and the window size (win_size) as the function arguments. As these 2 arguments are not in the form of array or list, the broadcast function will not iterate on them, but just pass them in each iteration.

function moving_avg(vec,win_size)
    n = size(vec)[1]
    return window_avg.(Ref(vec),win_size,1:n)[win_size:end]
end
function window_avg(a,win,i)
    if i>=win
        return sum(a[i-win+1:i])/win
    end
end

Regression using MLJ

Matthew leung — Tue, 07 Jun 2022 14:29:36 +0000

I try to analyze the daily price of Gold (GLD) with the relationship to the Daily Treasury Par Real Yield Curve Rates using the MLJ regression in Julia.

The Yield Curve Rates has 5-year rate, 7-year rate, 10-year rate, 20-year rate, and 30-year rate. I use all these rates, with its previous from 1 to 4 days values as the features to predict the price of Gold in the next day. Data from 2010–02–22 was collected, and split into training (70%) and testing data (30%). Here is the code in Julia/MLJ I used for the training and testing.

using MLJ
y, X = MLJ.unpack(full_df, ==(Symbol("target")), colname -> true)
train, test = partition(eachindex(y), 0.7)
using MLJLinearModels
reg = @load RidgeRegressor pkg = "MLJLinearModels"
model_reg = reg()
mach = machine(model_reg, X, y)
fit!(mach, rows=train)
yhat = predict(mach, X[test,:])

Ridge Regression Model is used for training and prediction. Here is the prediction value compared with the actual value.

Then I tried to use SHAP value to analysis the features importance.

using ShapML
function predict_function(model, data)
  data_pred = DataFrame(y_pred = predict(model, data))
  return data_pred
end
explain = copy(full_df[train, :]) 
# Remove the outcome column.
explain = select(explain, Not(Symbol("target"))) 
# An optional reference population to compute the baseline prediction. 
reference = copy(full_df)  
reference = select(reference, Not(Symbol("target")))
sample_size = 60  # Number of Monte Carlo samples.
data_shap = ShapML.shap(explain = explain,
                         reference = reference,
                         model = mach,
                         predict_function = predict_function,
                         sample_size = sample_size,
                         seed = 1
                         )

show(data_shap, allcols = true)

Then, plot the Mean Absolute Shapley Value with the Gadfly.jl

using Gadfly
data_plot = combine(DataFrames.groupby(data_shap, [:feature_name]),:shap_effect=>mean)
transform!(data_plot,:shap_effect_mean => x -> abs.(x))
data_plot = sort(data_plot, order(:shap_effect_mean_function, rev = true))

baseline = round(data_shap.intercept[1], digits = 1)

p = plot(data_plot, y = :feature_name, x = :shap_effect_mean_function, Coord.cartesian(yflip = true), Scale.y_discrete, Geom.bar(position = :dodge, orientation = :horizontal),Theme(bar_spacing = 1mm),Guide.xlabel("|Shapley effect| (baseline = $baseline)"), Guide.ylabel(nothing),Guide.title("Feature Importance - Mean Absolute Shapley Value"))
draw(PNG("features.png", 6inch, 6inch), p)

It can be seen that the 7-year yield rate is the most important features. Therefore, the Gold price is related to the 7-year yield rate. However, as the prediction is not fully matched with the actual value, there are other factors affecting the Gold price.

The above is originally posted on my blog in medium.com.

Use transformer to predict time series

Matthew leung — Thu, 02 Jun 2022 10:50:09 +0000

Just tried to use the transformer package to predict time series, details in https://iwasnothing.medium.com/use-transformer-in-julia-8409d2c0b642

Here is the summary:

//
The model can be coded easily using the 2 Julia packages: Flux, and Transformers.jl.

The full code can be found in my Github:
https://github.com/iwasnothing/julia_transformer_ts/blob/main/timeseries-transformer.jl

Data Preparation

Each training data sample is a sub-sequence of the time series by shifting 1 unit to the right.

Prediction Testing

For model prediction with the sub-sequence ix, the encoder input (ix[:,1:enc_seq_len,:]), and decoder input by shifting the target by 1 left (ix[:,enc_seq_len:sz[1]-1,:]) were fed into the encoder_forward and decoder_forward to predict the last value (dec[end,:]) in the sequence.

Test 1 with sin curve

First, I tried to verify the model with simple sin curve data, which it should be predicted well.

Both the predicted value and actual value are well fit each other. The result is good.

Test 2 with interest rate data

I would like to use the actual data to test the model. The data I use is the 10-year daily treasury real yield downloaded since 2003 from the US Treasury. I use 15-day moving average to smoothen the data.

Both the predicted value (blue) and actual value (red) are well fit each other. The result is good.

Test 3 with stock price data

Trying to predict the stock price which is more volatile than interest rate. I use the Julia Package MarketData to get the daily closing stock price. Again, I smoothen the data by taking 15-day moving average.

Both the predicted value (blue) and actual value (red) are not well fit each other. The result is not good.